analysis prey communities - prey biomass - size distribution Sticky traps were sampled in 2021 Each location was succesfull trap was present for 1 week
#data setup
datainput cleanup exploration
here()
## [1] "C:/Users/kadwolf/OneDrive - UGent/UGent-PC/Kadwolf/Documents/00_SPINCITY/henri/ST_mixed_model"
set.seed(4)
data.gent<-read_csv2("data_ST_regio_Gent.csv", na=c("NB",""))
data.antwerpen<-read_csv2("data_ST_regio_Antwerpen.csv", na=c("NB",""))
data.leuven<-read_csv2("data_ST_regio_Leuven.csv", na=c("NB",""))
#check if colnames are the same can be done via names() or colnames(df)==colnames(data2)
#colnames(df_gent)==colnames(df_antwerpen)
#colnames(df_gent)==colnames(df_leuven)
#colnames(df_antwerpen)==colnames(df_leuven)
#combine the data and put all the variables correctly
data.st<-bind_rows(data.gent, data.antwerpen, data.leuven)
dim(data.st)
## [1] 6386 20
remove extra sampling P20SE, and rename columns
### remove unneccessary columns: limval type person ST_subdiv and rename to english
data.st <- data.st %>%
rename("sampleID"="staalid",
"order"= "orde",
"suborder"="onderorde",
"length"="lengte",
"width"="breedte",
"remark"="opmerking") %>%
select(start_date, end_date, sampleID, project, plotid, U_landscape, U_local, location,
phylum, class, order, suborder, length, width, n_indiv, remark)
data.st <- data.st %>%
filter(location != "P20SE")
put all factors correctly
all 81 locations were sampled. P10SY and P13SY were less suitable (check if they are very different)
data.st$sampleID <- as.factor(data.st$sampleID)
data.st$project <- as.factor(data.st$project)
data.st$plotid <- as.factor(data.st$plotid)
levels(data.st$plotid) # 27 locations
## [1] "P01" "P02" "P03" "P04" "P05" "P06" "P07" "P08" "P09" "P10" "P11" "P12"
## [13] "P13" "P14" "P15" "P16" "P17" "P18" "P19" "P20" "P21" "P22" "P23" "P24"
## [25] "P25" "P26" "P27"
data.st$U_landscape <- as.factor(data.st$U_landscape)
data.st$U_landscape<- factor(data.st$U_landscape, levels=c("LOW", "MEDIUM", "HIGH"))
data.st$U_local <- as.factor(data.st$U_local)
data.st$U_local<- factor(data.st$U_local, levels=c("LOW", "MEDIUM", "HIGH"))
data.st$phylum <- factor(data.st$phylum)
data.st$class <- factor(data.st$class)
data.st$order <- factor(data.st$order)
data.st$suborder <- factor(data.st$suborder)
data.st$location <- factor(data.st$location)
data.st$location <- factor(data.st$location,c("P01SG", "P01SY", "P01SR", "P02SG", "P02SY", "P02SR", "P03SG", "P03SY", "P03SR", "P04SG", "P04SY", "P04SR", "P05SG", "P05SY", "P05SR", "P06SG", "P06SY", "P06SR", "P07SG", "P07SY", "P07SR", "P08SG", "P08SY", "P08SR", "P09SG", "P09SY", "P09SR", "P10SG", "P10SY", "P10SR", "P11SG", "P11SY", "P11SR", "P12SG", "P12SY", "P12SR", "P13SG", "P13SY", "P13SR", "P14SG", "P14SY", "P14SR", "P15SG", "P15SY", "P15SR", "P16SG", "P16SY", "P16SR", "P17SG", "P17SY", "P17SR", "P18SG", "P18SY", "P18SR", "P19SG", "P19SY", "P19SR", "P20SG", "P20SY", "P20SR", "P21SG", "P21SY", "P21SR", "P22SG", "P22SY", "P22SR", "P23SG", "P23SY", "P23SR", "P24SG", "P24SY", "P24SR", "P25SG", "P25SY", "P25SR", "P26SG", "P26SY", "P26SR", "P27SG", "P27SY", "P27SR"))
length(levels(data.st$location)) ## correct all 81 locations were successfully sampled
## [1] 81
data.st <- data.st %>%
mutate(urb_cat= paste0(U_landscape, U_local))
data.st$urb_cat <- as.factor(data.st$urb_cat)
data.st$urb_cat <- factor(data.st$urb_cat, levels=c("LOWLOW", "LOWMEDIUM", "LOWHIGH", "MEDIUMLOW", "MEDIUMMEDIUM", "MEDIUMHIGH", "HIGHLOW", "HIGHMEDIUM", "HIGHHIGH"))
urb.col9 <-c("#38A800", "lightgreen", "darkolivegreen3" , "yellow", "#E6E600", "goldenrod1", "deeppink3", "purple", "#8c008c") #colour vector for the 9 urbanisation categories
urb.col3<- c("#38A800", "#E6E600", "#8B008B") ##urbanisation colour vector to use in graphs
col_vec <- c(rep("#8B008B", 3), rep("#E6E600", 3), rep("#38A800", 3))
plotid.col27 <- rep(col_vec, 3) #handy to colour the facets boxes
data.st$day <- as.numeric(chron(dates. = as.character(data.st$start_date), format = c(dates = "d/m/y"), origin. = c(1, 1, 2021)))
## after this date can be put to a date object
data.st$start_date <- as.Date(data.st$start_date, format="%d/%m/%Y")
data.st$end_date <- as.Date(data.st$end_date, format="%d/%m/%Y")
glimpse(data.st)
## Rows: 6,282
## Columns: 18
## $ start_date <date> 2021-09-07, 2021-09-07, 2021-09-07, 2021-09-07, 2021-09-0…
## $ end_date <date> 2021-09-14, 2021-09-14, 2021-09-14, 2021-09-14, 2021-09-1…
## $ sampleID <fct> SC21P01SRLASEP07, SC21P01SYLASEP07, SC21P01SYLASEP07, SC21…
## $ project <fct> SC21, SC21, SC21, SC21, SC21, SC21, SC21, SC21, SC21, SC21…
## $ plotid <fct> P01, P01, P01, P01, P01, P01, P01, P01, P01, P01, P01, P01…
## $ U_landscape <fct> HIGH, HIGH, HIGH, HIGH, HIGH, HIGH, HIGH, HIGH, HIGH, HIGH…
## $ U_local <fct> HIGH, MEDIUM, MEDIUM, LOW, LOW, HIGH, MEDIUM, LOW, LOW, LO…
## $ location <fct> P01SR, P01SY, P01SY, P01SG, P01SG, P01SR, P01SY, P01SG, P0…
## $ phylum <fct> ARTHROPODA, ARTHROPODA, ARTHROPODA, ARTHROPODA, ARTHROPODA…
## $ class <fct> ARACHNIDA, ARACHNIDA, ARACHNIDA, INSECTA, INSECTA, INSECTA…
## $ order <fct> ARANEAE, ARANEAE, ARANEAE, COLEOPTERA, COLEOPTERA, COLEOPT…
## $ suborder <fct> NA, NA, NA, POLYPHAGA, POLYPHAGA, POLYPHAGA, POLYPHAGA, BR…
## $ length <dbl> 1.5, 2.0, 0.8, 1.1, 0.8, 2.6, 1.8, 0.5, 7.4, 7.5, 7.5, 6.4…
## $ width <dbl> 1.3, 0.9, 0.4, 0.5, 0.5, 0.7, 0.9, 0.1, 3.5, 3.1, 2.7, 2.7…
## $ n_indiv <dbl> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ remark <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ urb_cat <fct> HIGHHIGH, HIGHMEDIUM, HIGHMEDIUM, HIGHLOW, HIGHLOW, HIGHHI…
## $ day <dbl> 249, 249, 249, 249, 249, 249, 249, 249, 249, 249, 249, 249…
names(which(colSums(is.na(data.st)) > 0)) ## this is okay
## [1] "suborder" "remark"
for some location both sticky traps were measured, resolve this by dividing n_indiv
## for some locations both sticky traps were measured and the average should be taken of the n_indiv
data.st<- data.st %>%
mutate(sampleID = case_when(
sampleID == "SC21P05SGLASEP07" ~ "SC21P05SGLBSEP07",
sampleID == "SC21P06SGLAAUG28" ~ "SC21P06SGLBAUG28",
sampleID == "SC21P07SYLASEP07" ~ "SC21P07SYLBSEP07",
sampleID == "SC21P08SGLAAUG18" ~ "SC21P08SGLBAUG18",
sampleID == "SC21P09SGLAAUG28" ~ "SC21P09SGLBAUG28",
sampleID == "SC21P10SGLASEP08" ~ "SC21P10SGLBSEP08",
sampleID == "SC21P11SGLAAUG27" ~ "SC21P11SGLBAUG27",
sampleID == "SC21P12SYLAAUG17" ~ "SC21P12SYLBAUG17",
sampleID == "SC21P15SYLAAUG27" ~ "SC21P15SYLBAUG27",
sampleID == "SC21P17SRLAAUG17" ~ "SC21P17SRLBAUG17",
TRUE ~ sampleID
))
double_st <- c("SC21P05SGLBSEP07", "SC21P06SGLBAUG28", "SC21P07SYLBSEP07", "SC21P08SGLBAUG18", "SC21P09SGLBAUG28", "SC21P10SGLBSEP08", "SC21P11SGLBAUG27", "SC21P12SYLBAUG17", "SC21P15SYLBAUG27", "SC21P17SRLBAUG17")
data.st <- data.st %>%
mutate(n_indiv=(ifelse(sampleID %in% double_st, ceiling(n_indiv/2), n_indiv)))
next we will arrange the data in long so it is handy to use with ggplot
sum(data.st$n_indiv) #we want a dataset where there are as many rows as there are prey items measured, so equal to n_indiv
## [1] 9605
data.st <- data.st[rep(seq_len(nrow(data.st)), data.st$n_indiv), 1:ncol(data.st)]
nrow(data.st) # number of rows is now the same as the number of observations we only need to put the n_indiv to 1
## [1] 9605
data.st$n_indiv<-1
sum(data.st$n_indiv, na.rm=T) #double check if it is correct ok
## [1] 9605
calculate area of the insects the as measure for biomass
#create column with the surface area as a proxy of the biomass
#round the length measurements to 1 decimal
data.st$length <- round(data.st$length, 1)
data.st$width <- round(data.st$width, 1)
data.st <- data.st %>%
mutate(area=length*width)
#data exploration order
avg1 <- data.st %>%
group_by(plotid, U_landscape, U_local, phylum, class, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
avg2 <- data.st %>%
group_by(plotid, U_landscape, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
# for each urb land and urb local combination
avg3 <- data.st %>%
group_by(U_landscape, U_local, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
## U_local in columns
avg4 <-data.st %>%
group_by(U_landscape, U_local, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
pivot_wider(values_from = total_n_indiv, names_from=U_local)
## looks already like a kind of community matrix
avg5 <-data.st %>%
group_by(location, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
pivot_wider(values_from = total_n_indiv, names_from=location)
avg6 <-data.st %>%
group_by(U_landscape, U_local, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
pivot_wider(values_from = total_n_indiv, names_from=order)
sub_order
avg1_sub <- data.st %>%
group_by(plotid, U_landscape, U_local, phylum, class, order, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
avg2_sub <- data.st %>%
group_by(plotid, U_landscape, order, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
# for each urb land and urb local combination
avg3_sub <- data.st %>%
group_by(U_landscape, U_local, order, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T),
avg_length=mean(length,na.rm=T),
sd_avg_length = sd(length, na.rm=T))
## U_local in columns
avg4_sub <-data.st %>%
group_by(U_landscape, U_local, order, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
pivot_wider(values_from = total_n_indiv, names_from=U_local)
## looks already like a kind of community matrix
avg5_sub <-data.st %>%
group_by(location, order, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
pivot_wider(values_from = total_n_indiv, names_from=location)
data.st.tot <- data.st %>%
group_by(plotid, U_landscape, U_local, location) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))
fig_n_indiv <- ggplot(data.st.tot, aes(x = U_local, y = total_n_indiv, fill=U_local)) +
geom_col() +
scale_fill_manual(values=urb.col3, guide="none")+
facet_grid(~plotid)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens across each sampling location")+
theme_bw()+
theme(axis.text.x = element_text(angle=45,hjust=1, size=rel(0.8)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_n_indiv))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
data.st.tot[which.max(data.st.tot$total_n_indiv),]
## # A tibble: 1 × 5
## # Groups: plotid, U_landscape, U_local [1]
## plotid U_landscape U_local location total_n_indiv
## <fct> <fct> <fct> <fct> <dbl>
## 1 P09 LOW LOW P09SG 487
P09SG: largest number of specimens on the sticky trap
data.st.tot[which.min(data.st.tot$total_n_indiv),]
## # A tibble: 1 × 5
## # Groups: plotid, U_landscape, U_local [1]
## plotid U_landscape U_local location total_n_indiv
## <fct> <fct> <fct> <fct> <dbl>
## 1 P24 MEDIUM LOW P24SG 15
P24SG: smallest number of speciments on the sticky trap
##number of specimens along the urbanisation gradient
fig2 <- ggplot(data.st.tot, aes(x = U_local, y = total_n_indiv, color=U_local)) +
geom_point() +
scale_colour_manual(values=urb.col3, name="urbanisation level \nat local scale")+
facet_grid(~U_landscape)+
ylab("number of specimens")+
xlab("urbanisation at local scale")+
ggtitle("number of specimens along the urbanisation gradient")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1)),
plot.title = element_text(size=rel(1.1)))
fig2
g <- ggplot_gtable(ggplot_build(fig2))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
same graph: other visualisation :
ggplot(data.st.tot, aes(U_landscape, total_n_indiv, fill= U_local)) +
geom_boxplot(position=position_dodge(width=0.75))+
geom_jitter(position=position_dodge(width=0.75), alpha=0.6, size=1.5, color="gray3")+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
ylab("number of specimens")+
xlab("urbanisation level at landscape scale")+
ggtitle("number of specimens along the urbanisation gradient")+
theme_bw()+
theme(axis.text = element_text(size=rel(01)), axis.title=element_text(size=rel(1)))
explore order and suborder
##order ###by plotid graph showing the relative abundance of each arthropod order for each sampling location:
#16 ordes
col_ord <- c ("#323133", "#B8C0C7","#42A7F7", "blue4", "cyan", "goldenrod1", "#D9E33A", "yellow", "darkorchid4", "darkorchid1", "pink", "darkolivegreen1", "darkolivegreen4", "green", "firebrick1", "darkred" )
rel_n<- data.st %>%
group_by(plotid, U_landscape, U_local, location, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(plotid, U_landscape, U_local, location) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc ) %>%
ungroup()
### reorder the factors so it is plotted from highest number of specimens to lowest (more easy to compare)
rel_n$order <- reorder(rel_n$order, rel_n$rel_n_indiv)
rel_n$order <- factor(rel_n$order, levels=rev(levels(rel_n$order)))
## put also the order of data.st in this way (otherwise the colour code will not be preserved for each order
data.st$order <- factor(data.st$order, levels=levels(rel_n$order))
fig_rel_ord<- rel_n %>% ggplot(aes(x = U_local, y = rel_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~plotid)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of each taxonomic order for each sampling location")+
theme_bw()+
theme(axis.text.x = element_text(angle=45,hjust=1, size=rel(0.8)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_rel_ord))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
graph showing the number of specimens of each arthropod order for each sampling location:
ord_n <- data.st %>%
group_by(plotid, U_landscape, U_local, location, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()
fig_ord<- ord_n %>% ggplot(aes(x = U_local, y = total_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~plotid)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens of each taxonomic order for each sampling location")+
theme_bw()+
theme(axis.text.x = element_text(angle=45,hjust=1, size=rel(0.8)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_ord))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
the relative abundance of the Arthropod orders across different urbanisation levels at landscape and local scale:
rel_n_urb<- data.st %>%
group_by(U_landscape, U_local, urb_cat, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(U_landscape, U_local, urb_cat) %>%
mutate(n_indiv_urb = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_urb ) %>%
ungroup()
fig_rel_ord_urb<- rel_n_urb %>% ggplot(aes(x = U_local, y = rel_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~U_landscape)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of the taxonomic orders sampled along the urbanisation gradient")+
theme_bw()+
theme(axis.text.x = element_text(size=rel(1)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_rel_ord_urb))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
graph showing the number of specimens of each arthropod order across different urbanisation levels at landscape and local scale:
n_urb<- data.st %>%
group_by(U_landscape, U_local, urb_cat, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()
fig_n_ord_urb<- n_urb %>% ggplot(aes(x = U_local, y = total_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~U_landscape)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens of the taxonomic orders sampled along the urbanisation gradient")+
theme_bw()+
theme(axis.text.x = element_text(size=rel(1)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_n_ord_urb))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
##top3 orders diptera hemiptera hymenoptera are the top 3 orders
top3 <- c("DIPTERA", "HEMIPTERA", "HYMENOPTERA")
###by plotid graph showing the relative abundance of the top3 arthropod orders for each sampling location:
top3_rel_n<- data.st %>% filter(order %in% top3) %>%
group_by(plotid, U_landscape, U_local, location, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(plotid, U_landscape, U_local, location) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc ) %>%
ungroup()
fig_top3_rel_ord<- top3_rel_n %>% ggplot(aes(x = U_local, y = rel_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~plotid)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of the top3 athropod orders for each sampling location")+
theme_bw()+
theme(axis.text.x = element_text(angle=45,hjust=1, size=rel(0.8)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_top3_rel_ord))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
graph showing the number of specimens : graph showing the number of specimens of the top3 arthropod orders for each sampling location:
top3_ord_n <- data.st %>% filter(order %in% top3) %>%
group_by(plotid, U_landscape, U_local, location, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()
fig_top3_ord<- top3_ord_n %>% ggplot(aes(x = U_local, y = total_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~plotid)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens of the top3 arthropod orders for each location")+
theme_bw()+
theme(axis.text.x = element_text(angle=45,hjust=1, size=rel(0.8)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_top3_ord))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
###by urb level graph showing the relative abundance : the relative abundance of the top3 Arthropod orders across different urbanisation levels at landscape and local scale
top3_rel_n_urb<- data.st %>% filter(order %in% top3) %>%
group_by(U_landscape, U_local, urb_cat, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(U_landscape, U_local, urb_cat) %>%
mutate(n_indiv_urb = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_urb ) %>%
ungroup()
fig_top3_rel_ord_urb<- top3_rel_n_urb %>% ggplot(aes(x=U_local, y=rel_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~U_landscape)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of the top3 Arthropod orders sampled along the urbanisation gradient")+
theme_bw()+
theme(axis.text.x = element_text(size=rel(1)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_top3_rel_ord_urb))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
the same graph but with the number of specimens : the number of specimens of the top3 Arthropod orders across different urbanisation levels at landscape and local scale
top3_n_urb<- data.st %>% filter(order %in% top3) %>%
group_by(U_landscape, U_local, urb_cat, order) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()
fig_top3_n_ord_urb<- top3_n_urb %>% ggplot(aes(x = U_local, y = total_n_indiv, fill=order)) +
geom_col() +
scale_fill_manual(values=col_ord)+
facet_grid(~U_landscape)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens of the top3 Arthropod orders sampled along the urbanisation gradient")+
theme_bw()+
theme(axis.text.x = element_text(size=rel(1)),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(fig_top3_n_ord_urb))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
investigate suborder within the top3 orders
the number of specimens of the suborders within the top3 arthropod orders across different urbanisation levels at landscape and local scale
### for certain orders the suborder is not known so NA
## however for diptera, hemiptera and hymenoptera the data is complete
s1 <- data.st %>%
filter(order=="DIPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))%>%
ungroup()%>%
mutate(order="DIPTERA")
s2 <- data.st %>%
filter(order=="HEMIPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
mutate(order="HEMIPTERA")
s3 <- data.st %>%
filter(order=="HYMENOPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
mutate(order="HYMENOPTERA") %>%
ungroup()
s <- rbind(s1, s2, s3) ## number of specimens of the suborders within the top3 orders
fig_s<- s %>%
ggplot(aes(x = U_local, y = total_n_indiv, fill=suborder)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
facet_grid(order~plotid)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens per suborder of the top 3 orders for each location")+
theme_bw()+
theme(legend.position="bottom",
axis.text.x = element_text(angle=45,hjust=1, size=rel(0.7)),
plot.title = element_text(size=rel(1.1)))
col_ord_3 <- col_ord[1:3]
g1 <- ggplot_gtable(ggplot_build(fig_s))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(plotid.col27, 0.30), scales::alpha(col_ord_3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
the relative abundance of suborders within the top3 arthropod orders across different urbanisation levels at landscape and local scale:
r1<- data.st %>%
filter(order=="DIPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))%>%
ungroup()%>%
group_by(plotid, U_landscape, U_local) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="DIPTERA") %>%
ungroup()
r2<- data.st %>%
filter(order=="HEMIPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))%>%
ungroup()%>%
group_by(plotid, U_landscape, U_local) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="HEMIPTERA") %>%
ungroup()
r3<- data.st %>%
filter(order=="HYMENOPTERA") %>%
group_by(plotid, U_landscape, U_local, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))%>%
ungroup()%>%
group_by(plotid, U_landscape, U_local) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="HYMENOPTERA") %>%
ungroup()
r <- rbind(r1, r2, r3) ## relative abundance of the suborders within the top3 orders
fig_r <- r %>%
ggplot(aes(x = U_local, y = rel_n_indiv, fill=suborder)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
facet_grid(order~plotid)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of the suborders within the top 3 taxonomic orders for each sampling location")+
theme_bw()+
theme(legend.position="bottom",
axis.text.x = element_text(angle=45,hjust=1, size=rel(0.7)),
plot.title = element_text(size=rel(1.1)))
col_ord_3 <- col_ord[1:3]
g1 <- ggplot_gtable(ggplot_build(fig_r))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(plotid.col27, 0.30), scales::alpha(col_ord_3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
relative abundance of the suborders within the top3 orders for the different urbanisation levels at landscape and local scale:
### for certain orders the suborder is not known so NA
## however for diptera, hemiptera and hymenoptera the data is complete
rel_sub_1 <- data.st %>%
filter(order=="DIPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(U_landscape, U_local, urb_cat) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="DIPTERA") %>%
ungroup()
rel_sub_2 <- data.st %>%
filter(order=="HEMIPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(U_landscape, U_local, urb_cat) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="HEMIPTERA") %>%
ungroup()
rel_sub_3 <- data.st %>%
filter(order=="HYMENOPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
group_by(U_landscape, U_local, urb_cat) %>%
mutate(n_indiv_loc = sum(total_n_indiv, na.rm=T),
rel_n_indiv = total_n_indiv/n_indiv_loc,
order="HYMENOPTERA") %>%
ungroup()
rel_sub <- rbind(rel_sub_1, rel_sub_2, rel_sub_3)
fig_rel_subord<-rel_sub %>%
ggplot(aes(x = U_local, y = rel_n_indiv, fill=suborder)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
facet_grid(order~U_landscape)+
ylab("relative abundance")+
scale_y_continuous(labels=scales::percent_format())+
xlab("urbanisation level at local scale")+
ggtitle("relative abundance of the suborders within the top 3 orders \nacross different urbanisation levels at landscape and local scale")+
theme_bw()+
theme(plot.title = element_text(size=rel(1.1)))
col_ord_3 <- col_ord[1:3]
g1 <- ggplot_gtable(ggplot_build(fig_rel_subord))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(col_ord_3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
number of specimens of the suborders within the top3 orders for the different urbanisation levels at landscape and local scale
### for certain orders the suborder is not known so NA
## however for diptera, hemiptera and hymenoptera the data is complete
sub_1 <- data.st %>%
filter(order=="DIPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T))%>%
ungroup()%>%
mutate(order="DIPTERA")
sub_2 <- data.st %>%
filter(order=="HEMIPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
mutate(order="HEMIPTERA")
sub_3 <- data.st %>%
filter(order=="HYMENOPTERA") %>%
group_by(U_landscape, U_local, urb_cat, suborder) %>%
summarise(total_n_indiv=sum(n_indiv, na.rm=T)) %>%
ungroup()%>%
mutate(order="HYMENOPTERA") %>%
ungroup()
sub <- rbind(sub_1, sub_2, sub_3)
fig_subord<- sub %>%
ggplot(aes(x = U_local, y = total_n_indiv, fill=suborder)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
facet_grid(order~U_landscape)+
ylab("number of specimens")+
xlab("urbanisation level at local scale")+
ggtitle("number of specimens per suborder of the top 3 orders for each urbanisation level")+
theme_bw()+
theme(plot.title = element_text(size=rel(1.1)))
col_ord_3 <- col_ord[1:3]
g1 <- ggplot_gtable(ggplot_build(fig_subord))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(col_ord_3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
summary_top3 <- data.st%>%
filter(order %in% top3) %>%
droplevels()%>%
group_by(plotid, U_landscape, U_local, order, suborder) %>%
mutate(total_n_indiv=sum(n_indiv, na.rm=T))
### summary of the data with the total number of individuals per order)
#prey size distribution
##length
p<- ggplot(data.st, aes(x = U_local, y = length, color=U_local)) +
geom_point() +
scale_color_manual(values=urb.col3, name="urbanisation level at local scale")+
facet_grid(~plotid)+
ylab("length (mm)")+
theme(legend.position="bottom", axis.title.x=element_blank(),
axis.text.x = element_blank(), axis.ticks.x=element_blank(),
plot.title = element_text(size=rel(1.1)))
g <- ggplot_gtable(ggplot_build(p))
stript <- which(grepl('strip-t', g$layout$name))
fills <- scales::alpha(plotid.col27, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g$grobs[[i]]$grobs[[1]]$childrenOrder))
g$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g)
most individuals are small
## same histogram with ggplot
lenghist1 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=0.5) +
ggtitle("length per 0.5 mm") +
theme_bw()
lenghist2 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=0.25) +
ggtitle("length per 0.25 mm") +
theme_bw()
lenghist3 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=0.1)+
ggtitle("length per 0.10 mm") +
theme_bw()
ggarrange(lenghist1, lenghist2, lenghist3, ncol = 3)
area
areahist1 <- ggplot(data.st, aes(x=area)) +
geom_histogram(binwidth=1)+
ggtitle(bquote('area hist 1 ('~mm^2~')')) +
theme_bw()
areahist2 <- ggplot(data.st, aes(x=area)) +
geom_histogram(binwidth=0.1)+
ggtitle(bquote('area hist 0.1 ('~mm^2~')')) +
theme_bw()
areahist3 <- ggplot(data.st, aes(x=area)) +
geom_histogram(binwidth=0.05)+
ggtitle(bquote('area hist 0.05 ('~mm^2~')')) +
theme_bw()
ggarrange(areahist1, areahist2, areahist3, ncol = 3)
maybe best to use the length to divide the specimens in size classes? if we take all samples it looks like 0.1mm presents the best => but if we divide the samples in their urbanisation categories 0.5 or even 1 could work => then we have 45 size classes (otherwise 225)
hist1 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=0.1)+
ggtitle("histogram: length per 0.1 mm") +
facet_grid(U_landscape~U_local)+
scale_fill_manual(values=urb.col3)+
theme_bw()
hist1
g1 <- ggplot_gtable(ggplot_build(hist1))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(urb.col3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
hist2 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=0.5)+
ggtitle("histogram; length per 0.5 mm") +
facet_grid(U_landscape~U_local)+
scale_fill_manual(values=urb.col3)+
theme_bw()
hist2
g2 <- ggplot_gtable(ggplot_build(hist2))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(urb.col3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g2$grobs[[i]]$grobs[[1]]$childrenOrder))
g2$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g2)
hist3 <- ggplot(data.st, aes(x=length)) +
geom_histogram(binwidth=1)+
ggtitle("histogram: length per 1 mm") +
facet_grid(U_landscape~U_local)+
scale_fill_manual(values=urb.col3)+
theme_bw()
hist3
g3 <- ggplot_gtable(ggplot_build(hist3))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(urb.col3, 0.30))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g3$grobs[[i]]$grobs[[1]]$childrenOrder))
g3$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g3)
ggarrange(g1, g2, g3, nrow = 3)
##functional groups/species maybe work with 0.5 mm (or if necessary with 1mm but maybe this is to wide of a range). if we continue with 0.5 we have 45 size classes
breaks <-c(seq(0,22, by=0.5), Inf)
labels <- paste0("SC", 1:(length(breaks)-1))
labels <- c("SC01","SC02","SC03","SC04","SC05","SC06","SC07","SC08","SC09","SC10","SC11","SC12","SC13","SC14","SC15","SC16","SC17","SC18","SC19","SC20","SC21","SC22","SC23", "SC24", "SC25", "SC26", "SC27", "SC28", "SC29", "SC30", "SC31", "SC32", "SC33", "SC34", "SC35", "SC36", "SC37", "SC38", "SC39", "SC40", "SC41", "SC42", "SC43", "SC44", "SC45")
data.st$sizeclass <- cut(data.st$length, breaks = breaks, labels = labels, right = FALSE)
#then merge the data of the size classes with the taxonomy data
data.st$suborder <- ifelse(is.na(data.st$suborder), "ord", paste(data.st$suborder))
data.st <- data.st %>%
mutate(ord_sub = paste(str_sub(order, 1, 3), str_sub(suborder, 1, 3), sep="_"))
data.st <- data.st %>%
mutate(func_species = paste(sizeclass, ord_sub, sep="_"))
data.st$ord_sub <- as.factor(data.st$ord_sub)
length(levels(data.st$ord_sub))
## [1] 22
data.st$func_species <- as.factor(data.st$func_species)
length(levels(data.st$func_species)) # in total 176 functional species
## [1] 197
in total 197 functional species
sub_dip_col <- c("#31688EFF","#35B779FF")
sub_hem_col <- c("#443A83FF", "#21908CFF", "#8FD744FF")
sub_hym_col <- c("#440154FF", "#FDE725FF")
d<-data.st %>%
filter(order == "DIPTERA") %>%
droplevels() %>%
ggplot(aes(x = U_local, y = n_indiv, fill=func_species)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
theme(axis.text.x = element_text(angle=45,hjust=1)) +
ggtitle("DIPTERA functional groups") +
facet_grid(suborder~U_landscape)+
ylab("number of specimens per functional species")
g1 <- ggplot_gtable(ggplot_build(d))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(sub_dip_col, 0.50))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
d<-data.st %>%
filter(order == "HEMIPTERA") %>%
droplevels() %>%
ggplot(aes(x = U_local, y = n_indiv, fill=func_species)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
theme(axis.text.x = element_text(angle=45,hjust=1)) +
ggtitle("HEMIPTERA functional groups") +
facet_grid(suborder~U_landscape)+
ylab("number of specimens per functional species")
g1 <- ggplot_gtable(ggplot_build(d))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(sub_hem_col, 0.50))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
d<-data.st %>%
filter(order == "HYMENOPTERA") %>%
droplevels() %>%
ggplot(aes(x = U_local, y = n_indiv, fill=func_species)) +
geom_bar(stat="identity") +
scale_fill_viridis_d()+
theme(axis.text.x = element_text(angle=45,hjust=1)) +
ggtitle("HYMENOPTERA functional groups") +
facet_grid(suborder~U_landscape)+
ylab("number of specimens per functional species")
g1 <- ggplot_gtable(ggplot_build(d))
strip <- which(grepl('strip-', g1$layout$name))
fills <- c(scales::alpha(urb.col3, 0.30), scales::alpha(sub_hym_col, 0.50))
k <- 1
for (i in strip) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
dark blue or low size class number => small size individuals
#standardize the (continuous) explanatory variables so they have a mean of zero and standard deviation of one, so the estimated coefficients are all on the same scale and you can compare more easy the effect sizes
data.st$sday <- scale(data.st$day)
#data.st$slengte <- scale(data.st$length)
#data.st$sopp <- scale(data.st$area)
#assumptions mixed model: look at the residuals (they need to be +- normally distributed) plot(mixed.model, which=1 or 2) plot the residuals - the red line should be close to being flat, like the dashed grey line
poisson => but dharma test gives some strange results
no indication of correlation with sampling day U_local and U_landscape significant
df_N <- data.st %>%
group_by(plotid, U_landscape, U_local, urb_cat, sday, day) %>%
mutate(total_n_indiv=sum(n_indiv, na.rm=T)) %>% ungroup()
df <- df_N %>%
select(total_n_indiv, day) %>% ungroup()
df_cor<-cor(df, method = "pearson")
corrplot(df_cor,method = "circle", addgrid.col="grey", addCoef.col = "black", number.cex = 0.8, type="upper", title="correlation number of specimens and sampling day", mar=c(0,0,2,0))
df_N$U_landscape <- as.factor(df_N$U_landscape)
df_N$U_landscape<- factor(df_N$U_landscape, levels=c("LOW", "MEDIUM", "HIGH"))
df_N$U_local <- as.factor(df_N$U_local)
df_N$U_local<- factor(df_N$U_local, levels=c("LOW", "MEDIUM", "HIGH"))
modA<-glmmTMB(total_n_indiv ~ U_landscape*U_local+(1|plotid), data=df_N, family=poisson) #nbinom 2 to correct for overdispersion
summary(modA)
## Family: poisson ( log )
## Formula: total_n_indiv ~ U_landscape * U_local + (1 | plotid)
## Data: df_N
##
## AIC BIC logLik deviance df.resid
## 249584.8 249656.5 -124782.4 249564.8 9595
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plotid (Intercept) 0.1897 0.4356
## Number of obs: 9605, groups: plotid, 27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.139872 0.145211 35.40 <2e-16 ***
## U_landscapeMEDIUM -0.304649 0.205369 -1.48 0.138
## U_landscapeHIGH -0.140724 0.205364 -0.69 0.493
## U_localMEDIUM -0.228634 0.002851 -80.21 <2e-16 ***
## U_localHIGH -0.431410 0.003406 -126.68 <2e-16 ***
## U_landscapeMEDIUM:U_localMEDIUM 0.105811 0.004781 22.13 <2e-16 ***
## U_landscapeHIGH:U_localMEDIUM -0.107727 0.005177 -20.81 <2e-16 ***
## U_landscapeMEDIUM:U_localHIGH 0.301714 0.004892 61.67 <2e-16 ***
## U_landscapeHIGH:U_localHIGH 0.291675 0.005096 57.24 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3") #U_local and U_landscape is significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: total_n_indiv
## Chisq Df Pr(>Chisq)
## (Intercept) 1252.8705 1 <2e-16 ***
## U_landscape 2.2048 2 0.3321
## U_local 17325.0310 2 <2e-16 ***
## U_landscape:U_local 7005.2389 4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 2.0615, p-value = 0.064
## alternative hypothesis: two.sided
plot(allEffects(modA))
contrast(emmeans(modA, specs=~U_local), method="pairwise") |> as_tibble()
## # A tibble: 3 × 6
## contrast estimate SE df z.ratio p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LOW - MEDIUM 0.229 0.00215 Inf 107. 0
## 2 LOW - HIGH 0.234 0.00206 Inf 113. 0
## 3 MEDIUM - HIGH 0.00434 0.00237 Inf 1.83 0.160
modA_comp <- contrast(emmeans(modA, specs=~U_landscape*U_local), method="pairwise") |> as_tibble()
modA_comp %>% filter(p.value < 0.099) %>% knitr::kable()
| contrast | estimate | SE | df | z.ratio | p.value |
|---|---|---|---|---|---|
| LOW LOW - LOW MEDIUM | 0.2286342 | 0.0028505 | Inf | 80.20771 | 0 |
| LOW LOW - LOW HIGH | 0.4314100 | 0.0034055 | Inf | 126.67941 | 0 |
| MEDIUM LOW - MEDIUM MEDIUM | 0.1228229 | 0.0038389 | Inf | 31.99454 | 0 |
| MEDIUM LOW - MEDIUM HIGH | 0.1296961 | 0.0035122 | Inf | 36.92762 | 0 |
| HIGH LOW - HIGH MEDIUM | 0.3363613 | 0.0043218 | Inf | 77.82828 | 0 |
| HIGH LOW - HIGH HIGH | 0.1397353 | 0.0037906 | Inf | 36.86351 | 0 |
| LOW MEDIUM - LOW HIGH | 0.2027758 | 0.0035379 | Inf | 57.31556 | 0 |
| HIGH MEDIUM - HIGH HIGH | -0.1966259 | 0.0045611 | Inf | -43.10957 | 0 |
preds_modA <- emmeans(modA, specs=~U_landscape*U_local, type="response") |> as_tibble()
preds_modA
## # A tibble: 9 × 7
## U_landscape U_local rate SE df asymp.LCL asymp.UCL
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LOW LOW 171. 24.8 Inf 128. 227.
## 2 MEDIUM LOW 126. 18.3 Inf 94.7 167.
## 3 HIGH LOW 148. 21.5 Inf 112. 197.
## 4 LOW MEDIUM 136. 19.7 Inf 102. 181.
## 5 MEDIUM MEDIUM 111. 16.2 Inf 83.7 148.
## 6 HIGH MEDIUM 106. 15.4 Inf 79.7 141.
## 7 LOW HIGH 111. 16.1 Inf 83.4 147.
## 8 MEDIUM HIGH 111. 16.1 Inf 83.2 147.
## 9 HIGH HIGH 129. 18.7 Inf 97.0 171.
pd <- position_dodge(width=0.5)
g_N<-ggplot()+
geom_point(data=df_N,
aes(x=U_landscape, y=total_n_indiv,color=U_local),
position=pd, size=2)+
geom_pointrange(data=preds_modA,
aes(x=U_landscape, y=rate, ymin=asymp.LCL, ymax=asymp.UCL,
fill=U_local),
position=pd, size = 1, pch = 21)+
scale_fill_manual(values=urb.col3, name="Urbanisation level \nat local scale")+
scale_colour_manual(values=urb.col3, guide="none")+
labs(x="urbanisation level at landscape scale", y="number of specimens")+
ggtitle("effect of urbanisation on arthropod abundance")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)),
axis.title = element_text(size=rel(1.1)),
plot.title = element_text(size=rel(1.1)))
g_N
random factor : (1|plotid) fixed effects : urbanisation at landscape and local scale :U_landscape and U_local
sday : scaled sampling day test but did not influence the biomass data
calculated biomass as : W = a * L^b describing how body length (L) and biomass is correlated
a = 0.04 b = 2.26 source: Sabo et al. 2002: Sabo, J.L., Bastow, J.L. & Power, M.E., 2002. Length-Mass Relationships for Adult Aquatic and Terrestrial Invertebrates in a California Watershed.
(in study of the netherlands they us a very general a=0.03 b=2.63)
effect of urbanisation on arthropod biomass
data.st <- data.st %>%
mutate(biomass_indiv = 0.04*(length^2.26))
df_biomass <- data.st %>%
group_by(sampleID, plotid, U_landscape, U_local, urb_cat, day, sday) %>%
summarise(biomass = sum(biomass_indiv))%>% ungroup()
df <- df_biomass %>%
select(biomass, day) %>% ungroup()
df_cor<-cor(df, method = "pearson")
corrplot(df_cor,method = "circle", addgrid.col="grey", addCoef.col = "black", number.cex = 0.8, type="upper", title="correlation biomass and sampling day", mar=c(0,0,2,0))
modA<-glmmTMB(biomass ~ U_landscape+U_local+ (1|plotid),data=df_biomass, family=Gamma(link=log))
summary(modA)
## Family: Gamma ( log )
## Formula: biomass ~ U_landscape + U_local + (1 | plotid)
## Data: df_biomass
##
## AIC BIC logLik deviance df.resid
## 804.1 820.9 -395.1 790.1 74
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plotid (Intercept) 0.03642 0.1908
## Number of obs: 81, groups: plotid, 27
##
## Dispersion estimate for Gamma family (sigma^2): 0.42
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.4694 0.1744 25.627 < 2e-16 ***
## U_landscapeMEDIUM -0.2351 0.2083 -1.129 0.259066
## U_landscapeHIGH -0.3265 0.1997 -1.635 0.102153
## U_localMEDIUM -0.2088 0.1893 -1.103 0.270124
## U_localHIGH -0.6096 0.1841 -3.312 0.000927 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3") #U_local is significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: biomass
## Chisq Df Pr(>Chisq)
## (Intercept) 656.7630 1 < 2.2e-16 ***
## U_landscape 2.8112 2 0.245217
## U_local 11.5485 2 0.003107 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1348, p-value = 0.536
## alternative hypothesis: two.sided
plot(allEffects(modA))
contrast(emmeans(modA, specs=~U_local), method="pairwise") ## biomass different in low and high urbanised settings
## contrast estimate SE df z.ratio p.value
## LOW - MEDIUM 0.209 0.189 Inf 1.103 0.5123
## LOW - HIGH 0.610 0.184 Inf 3.312 0.0027
## MEDIUM - HIGH 0.401 0.181 Inf 2.209 0.0696
##
## Results are averaged over the levels of: U_landscape
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
preds_modA <- emmeans(modA, specs=~U_local, type="response") |> as_tibble()
preds_modA
## # A tibble: 3 × 6
## U_local response SE df asymp.LCL asymp.UCL
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LOW 72.4 9.55 Inf 55.9 93.8
## 2 MEDIUM 58.8 8.23 Inf 44.7 77.3
## 3 HIGH 39.4 5.31 Inf 30.2 51.3
g_modA<-ggplot() +
geom_pointrange(data=preds_modA,
aes(x=U_local, y=response, ymin=asymp.LCL, ymax=asymp.UCL, fill=U_local),
position=position_dodge(width=0.5), size = 1, pch = 21)+
scale_fill_manual(values=urb.col3, name="")+
labs(x="urbanisation level at local scale",
y="biomass", ylim=c(18,26))+
ggtitle("effect of urbanisation at local scale on arthropod biomass")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)),
axis.title = element_text(size=rel(1.1)),
plot.title = element_text(size=rel(1.1)))
g_modA
visualisation of statistics results of biomass:
g_modA_bis<-ggplot() +
geom_point(data=df_biomass, aes(x=U_local, y=biomass),
position=position_jitter(width = 0.2, height = 0), size=2, col="#B8C0C7")+
geom_pointrange(data=preds_modA,
aes(x=U_local, y=response, ymin=asymp.LCL, ymax=asymp.UCL, fill=U_local),
position=position_dodge(width=0.5), size = 1, pch = 21)+
scale_fill_manual(values=urb.col3, name="")+
labs(x="urbanisation level at local scale",
y=bquote('biomass'), ylim=c(18,26))+
ggtitle("effect of urbanisation level at local scale on arthropod biomass")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)),
axis.title = element_text(size=rel(1.1)),
plot.title = element_text(size=rel(1.1)))
g_modA_bis
other type of visualisation:
ggplot(df_biomass, aes(U_landscape, biomass, colour= U_local)) +
geom_point(position=position_dodge(width=0.40))+
scale_colour_manual(values = urb.col3, name="urbanisation \nat local scale") +
ylab(bquote('biomass '(mm^3)))+
xlab("urbanisation level at landscape scale")+
theme_bw()+
ggtitle("effect of urbanisation on arthropod biomass")+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1)))
other type of visualisation:
b<-ggplot(df_biomass, aes(U_landscape, biomass, fill= U_local)) +
geom_boxplot(position=position_dodge(width=0.75))+
geom_jitter(position=position_dodge(width=0.75), alpha=0.6, size=1.5, color="gray3")+
scale_fill_manual(values = urb.col3, name="urbanisation \nat local scale") +
ylab(bquote('biomass '(mm^3)))+
xlab("urbanisation level at landscape scale")+
theme_bw()+
ggtitle("effect of urbanisation on arthropod biomass")+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1)))
b
random factor : (1|plotid) fixed effects : urbanisation at landscape and local scale (interaction not significant)
response variable => log(length) to make it more normally distributed however levene test of residuals not ok
df <- data.st %>%
select(length, day) %>% ungroup()
df_cor<-cor(df, method = "pearson")
corrplot(df_cor,method = "circle", addgrid.col="grey", addCoef.col = "black", number.cex = 0.8, type="upper", title="correlation biomass and sampling period", mar=c(0,0,2,0))
#no correlation so we do not need to include day
modA<-glmmTMB(length ~ U_landscape+U_local+(1|plotid/location),data=data.st, family=gaussian)
summary(modA)
## Family: gaussian ( identity )
## Formula: length ~ U_landscape + U_local + (1 | plotid/location)
## Data: data.st
##
## AIC BIC logLik deviance df.resid
## 36496.6 36553.9 -18240.3 36480.6 9597
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## location:plotid (Intercept) 1.936e-01 0.4399602
## plotid (Intercept) 5.396e-07 0.0007346
## Residual 2.565e+00 1.6016842
## Number of obs: 9605, groups: location:plotid, 81; plotid, 27
##
## Dispersion estimate for gaussian family (sigma^2): 2.57
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.47491 0.11590 21.353 <2e-16 ***
## U_landscapeMEDIUM -0.06789 0.12884 -0.527 0.5982
## U_landscapeHIGH -0.06869 0.12814 -0.536 0.5919
## U_localMEDIUM 0.04068 0.12853 0.316 0.7516
## U_localHIGH -0.21957 0.12840 -1.710 0.0873 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: length
## Chisq Df Pr(>Chisq)
## (Intercept) 455.9548 1 <2e-16 ***
## U_landscape 0.3785 2 0.8276
## U_local 4.6931 2 0.0957 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9994, p-value = 0.992
## alternative hypothesis: two.sided
testZeroInflation(simulateResiduals(modA))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
plot(allEffects(modA))
contrast(emmeans(modA, specs=~U_local), method="pairwise")
## contrast estimate SE df t.ratio p.value
## LOW - MEDIUM -0.0407 0.129 9597 -0.316 0.9463
## LOW - HIGH 0.2196 0.128 9597 1.710 0.2014
## MEDIUM - HIGH 0.2602 0.129 9597 2.010 0.1099
##
## Results are averaged over the levels of: U_landscape
## P value adjustment: tukey method for comparing a family of 3 estimates
modB<-glmmTMB(log(length) ~ U_landscape+U_local+(1|plotid/location),data=data.st, family=gaussian)
summary(modB)
## Family: gaussian ( identity )
## Formula: log(length) ~ U_landscape + U_local + (1 | plotid/location)
## Data: data.st
##
## AIC BIC logLik deviance df.resid
## 16149.8 16207.2 -8066.9 16133.8 9597
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## location:plotid (Intercept) 2.106e-02 1.451e-01
## plotid (Intercept) 3.105e-09 5.572e-05
## Residual 3.087e-01 5.556e-01
## Number of obs: 9605, groups: location:plotid, 81; plotid, 27
##
## Dispersion estimate for gaussian family (sigma^2): 0.309
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.691153 0.038426 17.987 <2e-16 ***
## U_landscapeMEDIUM -0.020505 0.042781 -0.479 0.632
## U_landscapeHIGH 0.002699 0.042536 0.063 0.949
## U_localMEDIUM 0.001723 0.042664 0.040 0.968
## U_localHIGH -0.044774 0.042625 -1.050 0.294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modB, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log(length)
## Chisq Df Pr(>Chisq)
## (Intercept) 323.5202 1 <2e-16 ***
## U_landscape 0.3481 2 0.8403
## U_local 1.5111 2 0.4698
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modB))
testDispersion(simulateResiduals(modB))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0023, p-value = 0.872
## alternative hypothesis: two.sided
testZeroInflation(simulateResiduals(modB))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided
plot(allEffects(modB))
mod_test_gamma<-glmmTMB(length ~ U_landscape+U_local+(1|plotid/location),data=data.st, family=Gamma(link=log))
summary(mod_test_gamma)
## Family: Gamma ( log )
## Formula: length ~ U_landscape + U_local + (1 | plotid/location)
## Data: data.st
##
## AIC BIC logLik deviance df.resid
## 30179.2 30236.6 -15081.6 30163.2 9597
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## location:plotid (Intercept) 2.874e-02 0.1695297
## plotid (Intercept) 3.361e-08 0.0001833
## Number of obs: 9605, groups: location:plotid, 81; plotid, 27
##
## Dispersion estimate for Gamma family (sigma^2): 0.311
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.903554 0.044203 20.441 <2e-16 ***
## U_landscapeMEDIUM -0.040398 0.049035 -0.824 0.4100
## U_landscapeHIGH -0.033503 0.048802 -0.687 0.4924
## U_localMEDIUM 0.002405 0.048930 0.049 0.9608
## U_localHIGH -0.100438 0.048887 -2.054 0.0399 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_test_gamma, type="3")
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: length
## Chisq Df Pr(>Chisq)
## (Intercept) 417.8427 1 < 2e-16 ***
## U_landscape 0.7822 2 0.67631
## U_local 5.7096 2 0.05757 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(mod_test_gamma))
testDispersion(simulateResiduals(mod_test_gamma))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.3957, p-value < 2.2e-16
## alternative hypothesis: two.sided
testZeroInflation(simulateResiduals(mod_test_gamma))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
plot(allEffects(mod_test_gamma))
contrast(emmeans(mod_test_gamma, specs=~U_local), method="pairwise")
## contrast estimate SE df z.ratio p.value
## LOW - MEDIUM -0.00241 0.0489 Inf -0.049 0.9987
## LOW - HIGH 0.10044 0.0489 Inf 2.054 0.0995
## MEDIUM - HIGH 0.10284 0.0493 Inf 2.088 0.0923
##
## Results are averaged over the levels of: U_landscape
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
mod_test_pois<-glmmTMB(length ~ U_landscape+U_local+(1|plotid/location),data=data.st, family=poisson)
## Warning in glmmTMB(length ~ U_landscape + U_local + (1 | plotid/location), :
## non-integer counts in a poisson model
summary(mod_test_pois)
## Family: poisson ( log )
## Formula: length ~ U_landscape + U_local + (1 | plotid/location)
## Data: data.st
##
## AIC BIC logLik deviance df.resid
## 33215.3 33265.5 -16600.6 33201.3 9598
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## location:plotid (Intercept) 2.841e-02 0.1685532
## plotid (Intercept) 2.714e-08 0.0001648
## Number of obs: 9605, groups: location:plotid, 81; plotid, 27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.899999 0.044612 20.174 <2e-16 ***
## U_landscapeMEDIUM -0.037263 0.049806 -0.748 0.4544
## U_landscapeHIGH -0.031454 0.049473 -0.636 0.5249
## U_localMEDIUM 0.006068 0.049555 0.122 0.9025
## U_localHIGH -0.097941 0.049693 -1.971 0.0487 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_test_pois, type="3") #U_local is marginally significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: length
## Chisq Df Pr(>Chisq)
## (Intercept) 406.9794 1 < 2e-16 ***
## U_landscape 0.6550 2 0.72073
## U_local 5.4308 2 0.06618 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(mod_test_pois))
testDispersion(simulateResiduals(mod_test_pois))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.086, p-value = 0.024
## alternative hypothesis: two.sided
testZeroInflation(simulateResiduals(mod_test_pois))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
plot(allEffects(mod_test_pois))
mod_test_nbinom2<-glmmTMB(length ~ U_landscape+U_local+(1|plotid/location),data=data.st, family=poisson)
## Warning in glmmTMB(length ~ U_landscape + U_local + (1 | plotid/location), :
## non-integer counts in a poisson model
summary(mod_test_nbinom2)
## Family: poisson ( log )
## Formula: length ~ U_landscape + U_local + (1 | plotid/location)
## Data: data.st
##
## AIC BIC logLik deviance df.resid
## 33215.3 33265.5 -16600.6 33201.3 9598
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## location:plotid (Intercept) 2.841e-02 0.1685532
## plotid (Intercept) 2.714e-08 0.0001648
## Number of obs: 9605, groups: location:plotid, 81; plotid, 27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.899999 0.044612 20.174 <2e-16 ***
## U_landscapeMEDIUM -0.037263 0.049806 -0.748 0.4544
## U_landscapeHIGH -0.031454 0.049473 -0.636 0.5249
## U_localMEDIUM 0.006068 0.049555 0.122 0.9025
## U_localHIGH -0.097941 0.049693 -1.971 0.0487 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod_test_nbinom2, type="3") #U_local is marginally significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: length
## Chisq Df Pr(>Chisq)
## (Intercept) 406.9794 1 < 2e-16 ***
## U_landscape 0.6550 2 0.72073
## U_local 5.4308 2 0.06618 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(mod_test_nbinom2))
testDispersion(simulateResiduals(mod_test_nbinom2))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.086, p-value = 0.024
## alternative hypothesis: two.sided
testZeroInflation(simulateResiduals(mod_test_nbinom2))
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0, p-value < 2.2e-16
## alternative hypothesis: two.sided
plot(allEffects(mod_test_nbinom2))
statistics not yet fine: but gives that U_local is marginally significant (p = 0.06618) is probably the case => see visualisation :
### plot the distribution of the prey items
l1 <- data.st %>%
group_by(plotid, U_landscape) %>%
summarise(avg_length=mean(length, na.rm=T)) %>%
ggplot(aes(U_landscape, avg_length, fill= U_landscape)) +
geom_boxplot()+ geom_jitter(width=0.2)+
scale_fill_manual(values = urb.col3, guide="none") +
ylim(1.5,5)+
ylab("")+
ggtitle("landscape scale (3 km x 3 km)")+
theme_bw()+
theme(axis.title=element_text(size=rel(1.0)),
axis.text.x = element_text(size=rel(1.1)),
axis.text.y =element_text(size=rel(1.0)),
axis.title.x = element_blank())
l1
l2 <- data.st %>%
group_by(plotid, U_local) %>%
summarise(avg_length=mean(length, na.rm=T)) %>%
ggplot(aes(U_local, avg_length, fill= U_local)) +
geom_boxplot()+ geom_jitter(width=0.2)+
scale_fill_manual(values = urb.col3, guide="none")+
ylim(1.5,5)+
ylab("Average prey length (mm)")+
ggtitle("local scale (200 m x 200 m)")+
theme_bw()+
theme(axis.title.y=element_text(size=rel(1.1)),
axis.text.x =element_text(size=rel(1.1)),
axis.text.y =element_text(size=rel(1)),
axis.title.x=element_blank())
l2
fig1<- ggarrange(l2, l1, ncol = 2)
fig1
annotate_figure(fig1, bottom=text_grob("Urbanisation level", vjust=0.5, size=14),
top=text_grob("Effect of urbanisation on average arthropod prey size", vjust=0.2, size=15))
### plot the distribution of the prey items
l1 <- data.st %>%
group_by(plotid, U_landscape) %>%
summarise(avg_n_prey=mean(sum(n_indiv))) %>%
ggplot(aes(U_landscape, avg_n_prey, fill= U_landscape)) +
geom_boxplot()+
geom_jitter(width=0.2)+
scale_fill_manual(values = urb.col3, guide="none") +
xlab("")+
ylab("")+
ylim(100,800)+
ggtitle("landscape scale (3 km x 3 km)")+
theme_bw()+
theme(axis.title=element_text(size=rel(1.2)),
axis.text.x = element_text(size=rel(1.2)),
axis.text.y =element_text(size=rel(1.2)),
axis.title.x = element_blank(),
axis.title.y= element_blank())
l1
l2 <- data.st %>%
group_by(plotid, U_local) %>%
summarise(avg_n_prey=mean(sum(n_indiv))) %>%
ggplot(aes(U_local, avg_n_prey, fill= U_local)) +
geom_boxplot()+
geom_jitter(width=0.2)+
scale_fill_manual(values = urb.col3, guide="none") +
ylab("number of prey")+
ggtitle("local scale (200 m x 200 m)")+
theme_bw()+
theme(axis.title=element_text(size=rel(1.2)),
axis.text.x = element_text(size=rel(1.2)),
axis.text.y =element_text(size=rel(1.2)),
axis.title.x = element_blank())
l2
fig1<- ggarrange(l2, l1, ncol = 2)
fig1
annotate_figure(fig1, bottom=text_grob("Urbanization level", vjust=0.5, size=15))
all functional species
func.matrix <- data.st %>%
group_by(location, plotid, U_landscape, U_local, urb_cat, day, sday, func_species) %>%
summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat, day, sday)) %>%
column_to_rownames(var="location")
##dit is de community matrix
info <- data.st %>%
group_by(location, plotid, U_landscape, U_local, urb_cat, day, sday, func_species) %>%
summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat, day, sday)
## dit is de info
info
## # A tibble: 81 × 7
## # Groups: location, plotid, U_landscape, U_local, urb_cat, day, sday [81]
## location plotid U_landscape U_local urb_cat day sday[,1]
## <fct> <fct> <fct> <fct> <fct> <dbl> <dbl>
## 1 P01SG P01 HIGH LOW HIGHLOW 249 1.21
## 2 P01SY P01 HIGH MEDIUM HIGHMEDIUM 249 1.21
## 3 P01SR P01 HIGH HIGH HIGHHIGH 249 1.21
## 4 P02SG P02 HIGH LOW HIGHLOW 239 0.0798
## 5 P02SY P02 HIGH MEDIUM HIGHMEDIUM 239 0.0798
## 6 P02SR P02 HIGH HIGH HIGHHIGH 239 0.0798
## 7 P03SG P03 HIGH LOW HIGHLOW 229 -1.05
## 8 P03SY P03 HIGH MEDIUM HIGHMEDIUM 229 -1.05
## 9 P03SR P03 HIGH HIGH HIGHHIGH 229 -1.05
## 10 P04SG P04 MEDIUM LOW MEDIUMLOW 229 -1.05
## # ℹ 71 more rows
head(info)
## # A tibble: 6 × 7
## # Groups: location, plotid, U_landscape, U_local, urb_cat, day, sday [6]
## location plotid U_landscape U_local urb_cat day sday[,1]
## <fct> <fct> <fct> <fct> <fct> <dbl> <dbl>
## 1 P01SG P01 HIGH LOW HIGHLOW 249 1.21
## 2 P01SY P01 HIGH MEDIUM HIGHMEDIUM 249 1.21
## 3 P01SR P01 HIGH HIGH HIGHHIGH 249 1.21
## 4 P02SG P02 HIGH LOW HIGHLOW 239 0.0798
## 5 P02SY P02 HIGH MEDIUM HIGHMEDIUM 239 0.0798
## 6 P02SR P02 HIGH HIGH HIGHHIGH 239 0.0798
nmds_matrix=metaMDS(func.matrix, distance="bray", k=2, trymax=100, autotransform = F)
## Run 0 stress 0.205799
## Run 1 stress 0.2068969
## Run 2 stress 0.21016
## Run 3 stress 0.2127645
## Run 4 stress 0.2145304
## Run 5 stress 0.2071227
## Run 6 stress 0.2254975
## Run 7 stress 0.2035478
## ... New best solution
## ... Procrustes: rmse 0.02687038 max resid 0.1680316
## Run 8 stress 0.2089785
## Run 9 stress 0.210042
## Run 10 stress 0.2088447
## Run 11 stress 0.2052119
## Run 12 stress 0.2104413
## Run 13 stress 0.2051957
## Run 14 stress 0.2100039
## Run 15 stress 0.2059824
## Run 16 stress 0.2028529
## ... New best solution
## ... Procrustes: rmse 0.01817394 max resid 0.1503437
## Run 17 stress 0.2478003
## Run 18 stress 0.20595
## Run 19 stress 0.2101974
## Run 20 stress 0.2354425
## Run 21 stress 0.2100393
## Run 22 stress 0.2082859
## Run 23 stress 0.211463
## Run 24 stress 0.2085858
## Run 25 stress 0.2068969
## Run 26 stress 0.2076207
## Run 27 stress 0.2102914
## Run 28 stress 0.2234853
## Run 29 stress 0.2111056
## Run 30 stress 0.2202592
## Run 31 stress 0.2151924
## Run 32 stress 0.234025
## Run 33 stress 0.2402233
## Run 34 stress 0.2368402
## Run 35 stress 0.2082478
## Run 36 stress 0.2204409
## Run 37 stress 0.2326045
## Run 38 stress 0.2185666
## Run 39 stress 0.2278931
## Run 40 stress 0.2028529
## ... New best solution
## ... Procrustes: rmse 0.000122142 max resid 0.0009709275
## ... Similar to previous best
## *** Best solution repeated 1 times
nmds_matrix
##
## Call:
## metaMDS(comm = func.matrix, distance = "bray", k = 2, trymax = 100, autotransform = F)
##
## global Multidimensional Scaling using monoMDS
##
## Data: func.matrix
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2028529
## Stress type 1, weak ties
## Best solution was repeated 1 time in 40 tries
## The best solution was from try 40 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'func.matrix'
sppr<- specnumber(func.matrix) #species richness
sppr
## P01SG P01SY P01SR P02SG P02SY P02SR P03SG P03SY P03SR P04SG P04SY P04SR P05SG
## 29 28 28 38 33 20 24 35 31 25 33 31 57
## P05SY P05SR P06SG P06SY P06SR P07SG P07SY P07SR P08SG P08SY P08SR P09SG P09SY
## 33 25 45 26 23 46 42 28 49 36 22 56 36
## P09SR P10SG P10SY P10SR P11SG P11SY P11SR P12SG P12SY P12SR P13SG P13SY P13SR
## 29 53 20 29 49 31 28 40 46 26 37 11 27
## P14SG P14SY P14SR P15SG P15SY P15SR P16SG P16SY P16SR P17SG P17SY P17SR P18SG
## 43 39 41 30 54 24 26 34 32 37 44 42 44
## P18SY P18SR P19SG P19SY P19SR P20SG P20SY P20SR P21SG P21SY P21SR P22SG P22SY
## 27 28 32 26 20 16 28 30 26 26 36 23 22
## P22SR P23SG P23SY P23SR P24SG P24SY P24SR P25SG P25SY P25SR P26SG P26SY P26SR
## 32 31 33 19 13 17 21 27 23 26 35 36 34
## P27SG P27SY P27SR
## 26 23 28
df_species_rich <- sppr %>%
enframe() %>%
full_join(info, by = c("name"="location")) %>%
rename("species_richness" = "value")
ggplot(df_species_rich, aes(x = U_landscape, y = species_richness, fill = U_local)) +
geom_boxplot(position=position_dodge(width=0.75)) +
geom_point(position=position_dodge(width=0.75), colour="gray3")+
scale_colour_manual(values = urb.col3) +
scale_fill_manual(values=urb.col3)+
theme_bw()+
#theme(legend.position = "none",
labs(x = "urbanisation level at local scale",
y = "Number of functional species per subplot",
title = "Species richness")
number of species decreases with increasing urbanisation at local
scale
###statistics species richness species richness : started with poisson but overdispersion so modified to nbinm2 => U_local significant
df <- df_species_rich %>%
select(species_richness, day) %>% ungroup()
df_cor<-cor(df, method = "pearson")
corrplot(df_cor,method = "circle", addgrid.col="grey", addCoef.col = "black", number.cex = 0.8, type="upper", title="correlation species richness and sampling day", mar=c(0,0,2,0))
modA<-glmmTMB(species_richness ~ U_landscape+U_local+(1|plotid), data=df_species_rich, family=nbinom2) #nbinom 2 to correct for overdispersion
summary(modA)
## Family: nbinom2 ( log )
## Formula: species_richness ~ U_landscape + U_local + (1 | plotid)
## Data: df_species_rich
##
## AIC BIC logLik deviance df.resid
## 594.0 610.8 -290.0 580.0 74
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plotid (Intercept) 0.00804 0.08966
## Number of obs: 81, groups: plotid, 27
##
## Dispersion parameter for nbinom2 family (): 25
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.63379 0.07171 50.67 < 2e-16 ***
## U_landscapeMEDIUM -0.11987 0.08435 -1.42 0.15528
## U_landscapeHIGH -0.09885 0.08411 -1.18 0.23991
## U_localMEDIUM -0.12711 0.07222 -1.76 0.07839 .
## U_localHIGH -0.22610 0.07320 -3.09 0.00201 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3") #U_local is significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: species_richness
## Chisq Df Pr(>Chisq)
## (Intercept) 2567.5276 1 < 2e-16 ***
## U_landscape 2.3217 2 0.31323
## U_local 9.6391 2 0.00807 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0356, p-value = 0.768
## alternative hypothesis: two.sided
plot(allEffects(modA))
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
contrast(emmeans(modA, specs=~U_local), method="pairwise")
## contrast estimate SE df z.ratio p.value
## LOW - MEDIUM 0.127 0.0722 Inf 1.760 0.1832
## LOW - HIGH 0.226 0.0732 Inf 3.089 0.0057
## MEDIUM - HIGH 0.099 0.0741 Inf 1.336 0.3754
##
## Results are averaged over the levels of: U_landscape
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
preds_modA <- emmeans(modA, specs=~U_local, type="response") |> as_tibble()
preds_modA
## # A tibble: 3 × 6
## U_local response SE df asymp.LCL asymp.UCL
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 LOW 35.2 1.89 Inf 31.7 39.1
## 2 MEDIUM 31.0 1.70 Inf 27.8 34.5
## 3 HIGH 28.1 1.56 Inf 25.2 31.3
g_sp<-ggplot()+
geom_point(data=df_species_rich, aes(x=U_local, y=species_richness),
position=position_jitter(width = 0.2, height = 0), size=2, col="#B8C0C7")+
geom_pointrange(data=preds_modA,
aes(x=U_local, y=response, ymin=asymp.LCL, ymax=asymp.UCL, fill=U_local), position=position_dodge(width=0.5), size = 1, pch = 21)+
scale_fill_manual(values=urb.col3, name="")+
labs(x="urbanisation level at local scale", y="species richness")+
ggtitle("effect of urbanisation on species richness")+
theme_bw()
g_sp
NMDS
## new way:
NMDS.vect <- scores(nmds_matrix, display = "species") %>%
as.data.frame()
NMDS.scores <- scores(nmds_matrix, display = "sites") %>%
as.data.frame()%>%
rownames_to_column("location") %>%
full_join(info, by="location")
stressplot(nmds_matrix)
##define a theme to use later
clean_background <- theme(plot.background = element_rect("white"),
panel.background = element_rect("white"),
panel.grid = element_line("white"),
axis.line = element_line("gray25"),
axis.text = element_text(size = 12, color = "gray25"),
axis.title = element_text(color = "gray25"),
legend.text = element_text(size = 12),
legend.key = element_rect("white"))
ordbray.scores <- as.data.frame(nmds_matrix$points) %>%
bind_cols(info)
ordbray.scores$U_local <- factor(ordbray.scores$U_local, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.scores$U_landscape <- factor(ordbray.scores$U_landscape, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.hulls <- ordbray.scores %>%
group_by(U_local) %>%
slice(chull(MDS1, MDS2))
# plot nMDS with color and convex hulls as U_local
ordbray.plot.mds <- ggplot() +
geom_point(data = ordbray.scores, aes(x = MDS1, y = MDS2, color = U_local), size = 2) +
geom_polygon(data = ordbray.hulls, aes(x = MDS1, MDS2, fill = U_local), alpha = 0.5) +
scale_fill_manual(values=urb.col3) +
scale_colour_manual(values=urb.col3)+
theme_bw() +
xlab("NMDS1") +
ylab("NMDS2") +
ggtitle("ordination functional species : nmds bray curtis dissimilarity ")
ordbray.plot.mds
envfit on func species
# envfit() takes the output of metaMDS() and the species matrix you created
fit <- envfit(nmds_matrix, func.matrix, perm = 999)
#extract p-values for each species
fit_pvals <- fit$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_spp <- fit %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
full_join(., fit_pvals, by = "func_species") %>%
filter(pvals <= 0.001)
set.seed(4)
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_local, shape = U_landscape), size = 3, alpha = 0.8) +
stat_ellipse(aes(color = U_local)) +
scale_colour_manual(name= "Urbanisation \nat local scale", values= urb.col3)+
scale_shape_manual(name= "urbanisation level \nat landscape scale", values=c(16, 17, 15))+
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.3, height=0.1)) +
clean_background
nmds_plot_new
##permanova full dataset
library(pairwiseAdonis)
ST.dist<-vegdist(func.matrix, distance="bray")
perman <- adonis2(ST.dist~ plotid + U_landscape * U_local,
data=info,
strata=info$plotid,
permutations=999)
perman
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = ST.dist ~ plotid + U_landscape * U_local, data = info, permutations = 999, strata = info$plotid)
## Df SumOfSqs R2 F Pr(>F)
## plotid 26 7.3473 0.40003 1.4122 0.09 .
## U_local 2 0.6025 0.03280 1.5054 0.03 *
## U_landscape:U_local 4 0.8118 0.04420 1.0142 0.44
## Residual 48 9.6052 0.52297
## Total 80 18.3668 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(perman)
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| plotid | 26 | 7.3472568 | 0.4000300 | 1.412166 | 0.09 |
| U_local | 2 | 0.6024680 | 0.0328021 | 1.505352 | 0.03 |
| U_landscape:U_local | 4 | 0.8118197 | 0.0442005 | 1.014223 | 0.44 |
| Residual | 48 | 9.6052182 | 0.5229674 | NA | NA |
| Total | 80 | 18.3667627 | 1.0000000 | NA | NA |
bd_local <- betadisper(ST.dist, info$U_local)
bd_local
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = ST.dist, group = info$U_local)
##
## No. of Positive Eigenvalues: 58
## No. of Negative Eigenvalues: 22
##
## Average distance to median:
## LOW MEDIUM HIGH
## 0.4536 0.4498 0.4743
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 80 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 2.8477 2.0489 1.1139 0.8832 0.8674 0.8456 0.6520 0.6361
print(permutest(bd_local,pairwise=TRUE,permutations=999)) ## als het niet significant is, dan is het goed
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.00941 0.0047038 0.4998 999 0.615
## Residuals 78 0.73402 0.0094105
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.89100 0.460
## MEDIUM 0.88872 0.318
## HIGH 0.45456 0.32280
bd_land <- betadisper(ST.dist, info$U_landscape)
print(permutest(bd_land,pairwise=TRUE,permutations=999))
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.02302 0.0115095 1.4515 999 0.226
## Residuals 78 0.61848 0.0079292
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.242000 0.070
## MEDIUM 0.254571 0.675
## HIGH 0.075275 0.691879
no signs of “dispersion” => betadisper test is ok
info$interaction_comb <- interaction(info$U_landscape, info$U_local)
info$interaction_comb <- as.factor(info$interaction_comb)
pairwise.adonis(ST.dist, info$U_local) |> as_tibble() |> knitr::kable()
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| LOW vs MEDIUM | 1 | 0.3389884 | 1.5344068 | 0.0286621 | 0.069 | 0.207 | |
| LOW vs HIGH | 1 | 0.3448402 | 1.4804737 | 0.0276825 | 0.085 | 0.255 | |
| MEDIUM vs HIGH | 1 | 0.2198734 | 0.9585059 | 0.0180992 | 0.509 | 1.000 |
pairwise.adonis(ST.dist, info$U_landscape)
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 HIGH vs MEDIUM 1 0.1638342 0.6817194 0.01294034 0.882 1.000
## 2 HIGH vs LOW 1 0.3235652 1.4408247 0.02696112 0.078 0.234
## 3 MEDIUM vs LOW 1 0.2672431 1.2028779 0.02260926 0.197 0.591
no significant difference in communities
data.st.dhh <- data.st %>%
dplyr::filter(order %in% top3 & location !="P24SG" & location !="P20SG") %>%
droplevels()
dhh.matrix <- data.st.dhh %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat)) %>%
column_to_rownames(var="location")
##dit is de community matrix
##sum(func.matrix) gives now 8812
info <- data.st.dhh %>%
dplyr::group_by(location, plotid, U_landscape, U_local,urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat)
## dit is de info
#func.hell <- decostand(dhh.matrix, method = "hellinger") #veranderde niets
nmds_matrix=metaMDS(dhh.matrix, distance="bray", k=2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2931411
## Run 1 stress 0.2948932
## Run 2 stress 0.3022437
## Run 3 stress 0.3072661
## Run 4 stress 0.296218
## Run 5 stress 0.3004433
## Run 6 stress 0.3043133
## Run 7 stress 0.3032487
## Run 8 stress 0.3058203
## Run 9 stress 0.3093865
## Run 10 stress 0.2934998
## ... Procrustes: rmse 0.0285168 max resid 0.1294497
## Run 11 stress 0.3072166
## Run 12 stress 0.2986092
## Run 13 stress 0.3067337
## Run 14 stress 0.296813
## Run 15 stress 0.3143007
## Run 16 stress 0.3055637
## Run 17 stress 0.3064587
## Run 18 stress 0.2978105
## Run 19 stress 0.2964639
## Run 20 stress 0.2977014
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 5: no. of iterations >= maxit
## 15: stress ratio > sratmax
nmds_matrix
##
## Call:
## metaMDS(comm = dhh.matrix, distance = "bray", k = 2)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dhh.matrix))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2931411
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 0 (metric scaling or null solution)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dhh.matrix))'
sppr<- specnumber(dhh.matrix)
sppr
## P01SG P01SY P01SR P02SG P02SY P02SR P03SG P03SY P03SR P04SG P04SY P04SR P05SG
## 21 19 23 30 25 19 16 32 25 22 26 28 45
## P05SY P05SR P06SG P06SY P06SR P07SG P07SY P07SR P08SG P08SY P08SR P09SG P09SY
## 27 20 40 22 17 37 37 24 37 27 18 44 28
## P09SR P10SG P10SY P10SR P11SG P11SY P11SR P12SG P12SY P12SR P13SG P13SY P13SR
## 22 44 17 21 40 27 24 34 37 23 31 11 24
## P14SG P14SY P14SR P15SG P15SY P15SR P16SG P16SY P16SR P17SG P17SY P17SR P18SG
## 36 31 34 26 43 22 24 31 29 31 34 34 36
## P18SY P18SR P19SG P19SY P19SR P20SY P20SR P21SG P21SY P21SR P22SG P22SY P22SR
## 23 25 29 20 17 25 28 21 20 32 20 20 26
## P23SG P23SY P23SR P24SY P24SR P25SG P25SY P25SR P26SG P26SY P26SR P27SG P27SY
## 28 27 15 14 16 25 20 18 29 25 28 23 18
## P27SR
## 26
df_species_rich <- sppr %>%
enframe() %>%
full_join(info, by = c("name"="location")) %>%
rename("species_richness" = "value")
#modA<-glmmTMB(species_richness ~ U_landscape+U_local+(1|plotid), data=df_species_rich, family=nbinom2) #nbinom 2 to correct for overdispersion
#summary(modA)
#Anova(modA, type="3") #U_local is significant
#plot(simulateResiduals(modA))
#testDispersion(simulateResiduals(modA))
#plot(allEffects(modA))
top3_rich<-ggplot(df_species_rich, aes(x = U_landscape, y = species_richness, fill = U_local)) +
geom_boxplot(position=position_dodge(width=0.75)) +
geom_point(position=position_dodge(width=0.75), colour="gray3")+
scale_colour_manual(values = urb.col3) +
scale_fill_manual(values=urb.col3)+
theme_bw()+
#theme(legend.position = "none",
labs(x = "urbanisation level at local scale",
y = "Number of functional species per subplot",
title = "Species richness funct species top 3 orders")
#top3_rich ## silence for now (give all the same trend as the full dataset)
## new way:
NMDS.vect <- scores(nmds_matrix, display = "species") %>%
as.data.frame()
NMDS.scores <- scores(nmds_matrix, display = "sites") %>%
as.data.frame()%>%
rownames_to_column("location") %>%
full_join(info, by="location")
stressplot(nmds_matrix)
# envfit() takes the output of metaMDS() and the species matrix you created
fit <- envfit(nmds_matrix, dhh.matrix, perm = 999)
#extract p-values for each species
fit_pvals <- fit$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_spp <- fit %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
full_join(., fit_pvals, by = "func_species") %>%
filter(pvals <= 0.001)
set.seed(4)
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_local, shape = U_landscape), size = 3, alpha = 0.8) +
stat_ellipse(aes(color = U_local)) +
scale_color_manual(name="urbanisation level \nat local scale", values = urb.col3) +
scale_shape_manual(name= "urbanisation level \nat landscape scale", values=c(16, 17, 15))+
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.3, height=0.1)) +
clean_background
nmds_plot_new
ordbray.scores <- as.data.frame(nmds_matrix$points) %>%
bind_cols(info)
ordbray.scores$U_local <- factor(ordbray.scores$U_local, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.scores$U_landscape <- factor(ordbray.scores$U_landscape, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.hulls <- ordbray.scores %>%
group_by(U_local) %>%
slice(chull(MDS1, MDS2))
# plot nMDS with color and convex hulls as U_local
ordbray.plot.mds <- ggplot() +
geom_point(data = ordbray.scores, aes(x = MDS1, y = MDS2, color = U_local), size = 2) +
geom_polygon(data = ordbray.hulls, aes(x = MDS1, MDS2, fill = U_local), alpha = 0.5) +
scale_fill_manual(values=urb.col3) +
scale_colour_manual(values=urb.col3)+
theme_bw() +
xlab("NMDS1") +
ylab("NMDS2") +
ggtitle("ordination top3 order functional species : nmds bray curtis dissimilarity ")
ordbray.plot.mds
ordbray.scores <- as.data.frame(nmds_matrix$points) %>%
bind_cols(info)
ordbray.scores$U_local <- factor(ordbray.scores$U_local, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.scores$U_landscape <- factor(ordbray.scores$U_landscape, levels=c("LOW", "MEDIUM", "HIGH"))
ordbray.hulls <- ordbray.scores %>%
group_by(urb_cat) %>%
slice(chull(MDS1, MDS2))
# plot nMDS with color and convex hulls as U_local
ordbray.plot.mds <- ggplot() +
geom_point(data = ordbray.scores, aes(x = MDS1, y = MDS2, color = urb_cat), size = 2) +
geom_polygon(data = ordbray.hulls, aes(x = MDS1, MDS2, fill = urb_cat), alpha = 0.5) +
scale_fill_manual(values=urb.col9) +
scale_colour_manual(values=urb.col9)+
theme_bw() +
xlab("NMDS1") +
ylab("NMDS2") +
ggtitle("ordination top3 order functional species : nmds bray curtis dissimilarity ")
ordbray.plot.mds
library(pairwiseAdonis)
ST.dist<-vegdist(dhh.matrix, distance="bray")
perman <- adonis2(ST.dist~ plotid + U_landscape * U_local,
data=info,
strata=info$plotid,
permutations=999)
perman
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = ST.dist ~ plotid + U_landscape * U_local, data = info, permutations = 999, strata = info$plotid)
## Df SumOfSqs R2 F Pr(>F)
## plotid 26 6.8955 0.41449 1.4555 0.054 .
## U_local 2 0.6517 0.03918 1.7885 0.008 **
## U_landscape:U_local 4 0.7075 0.04253 0.9707 0.571
## Residual 46 8.3815 0.50381
## Total 78 16.6362 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(perman)
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| plotid | 26 | 6.8954858 | 0.4144864 | 1.4555461 | 0.054 |
| U_local | 2 | 0.6517387 | 0.0391759 | 1.7884547 | 0.008 |
| U_landscape:U_local | 4 | 0.7074600 | 0.0425253 | 0.9706806 | 0.571 |
| Residual | 46 | 8.3815314 | 0.5038124 | NA | NA |
| Total | 78 | 16.6362158 | 1.0000000 | NA | NA |
bd_local <- betadisper(ST.dist, info$U_local)
bd_local
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = ST.dist, group = info$U_local)
##
## No. of Positive Eigenvalues: 54
## No. of Negative Eigenvalues: 24
##
## Average distance to median:
## LOW MEDIUM HIGH
## 0.4172 0.4391 0.4631
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 78 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 2.9265 2.0659 1.1090 0.9075 0.8422 0.7317 0.6595 0.6253
print(permutest(bd_local,pairwise=TRUE,permutations=999)) ## als het niet significant is, dan is het goed
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.02748 0.0137417 1.5098 999 0.229
## Residuals 76 0.69172 0.0091016
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.417000 0.098
## MEDIUM 0.400740 0.365
## HIGH 0.094216 0.361372
bd_land <- betadisper(ST.dist, info$U_landscape)
print(permutest(bd_land,pairwise=TRUE,permutations=999))
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.01743 0.0087134 1.0954 999 0.37
## Residuals 76 0.60452 0.0079542
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.54500 0.119
## MEDIUM 0.52813 0.460
## HIGH 0.12164 0.42801
no signs of “dispersion” => betadisper test is ok
info$interaction_comb <- interaction(info$U_landscape, info$U_local)
info$interaction_comb <- as.factor(info$interaction_comb)
pairwise.adonis(ST.dist, info$U_local) |> as_tibble() |> knitr::kable()
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| LOW vs MEDIUM | 1 | 0.3556410 | 1.7858194 | 0.0344847 | 0.047 | 0.141 | |
| LOW vs HIGH | 1 | 0.4084134 | 1.9312649 | 0.0371889 | 0.023 | 0.069 | |
| MEDIUM vs HIGH | 1 | 0.2166091 | 0.9846406 | 0.0185835 | 0.468 | 1.000 |
pairwise.adonis(ST.dist, info$U_landscape)
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 HIGH vs MEDIUM 1 0.1397605 0.6311697 0.01246603 0.901 1.000
## 2 HIGH vs LOW 1 0.2979443 1.4014190 0.02674391 0.119 0.357
## 3 MEDIUM vs LOW 1 0.2287097 1.1143197 0.02138222 0.305 0.915
data.st.dip <- data.st %>%
dplyr::filter(order == "DIPTERA") %>%
droplevels()
dip.matrix <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat)) %>%
column_to_rownames(var="location")
##dit is de community matrix
##sum(dip.matrix) gives now 5287
info <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat)
## dit is de info
nmds_matrix=metaMDS(dip.matrix, distance="bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2647948
## Run 1 stress 0.2802532
## Run 2 stress 0.2759063
## Run 3 stress 0.2692734
## Run 4 stress 0.2718641
## Run 5 stress 0.2776076
## Run 6 stress 0.2840558
## Run 7 stress 0.2647879
## ... New best solution
## ... Procrustes: rmse 0.003919021 max resid 0.02634658
## Run 8 stress 0.2814574
## Run 9 stress 0.2793722
## Run 10 stress 0.2698903
## Run 11 stress 0.2811248
## Run 12 stress 0.2693693
## Run 13 stress 0.2758939
## Run 14 stress 0.2822895
## Run 15 stress 0.2786681
## Run 16 stress 0.2687636
## Run 17 stress 0.2794938
## Run 18 stress 0.2766532
## Run 19 stress 0.2751782
## Run 20 stress 0.282888
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmax
nmds_matrix
##
## Call:
## metaMDS(comm = dip.matrix, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dip.matrix))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2647879
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 7 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dip.matrix))'
sppr<- specnumber(dip.matrix)
sppr
## P01SG P01SY P01SR P02SG P02SY P02SR P03SG P03SY P03SR P04SG P04SY P04SR P05SG
## 13 10 10 19 13 10 9 20 14 14 10 18 26
## P05SY P05SR P06SG P06SY P06SR P07SG P07SY P07SR P08SG P08SY P08SR P09SG P09SY
## 12 10 22 13 10 22 19 12 21 18 7 25 17
## P09SR P10SG P10SY P10SR P11SG P11SY P11SR P12SG P12SY P12SR P13SG P13SY P13SR
## 13 20 7 8 20 17 12 18 18 12 15 5 8
## P14SG P14SY P14SR P15SG P15SY P15SR P16SG P16SY P16SR P17SG P17SY P17SR P18SG
## 18 18 20 15 22 14 12 16 14 19 21 19 21
## P18SY P18SR P19SG P19SY P19SR P20SG P20SY P20SR P21SG P21SY P21SR P22SG P22SY
## 13 12 17 9 7 7 14 13 8 11 13 12 11
## P22SR P23SG P23SY P23SR P24SG P24SY P24SR P25SG P25SY P25SR P26SG P26SY P26SR
## 8 13 15 3 7 9 8 15 11 10 17 16 17
## P27SG P27SY P27SR
## 13 9 14
df_species_rich <- sppr %>%
enframe() %>%
full_join(info, by = c("name"="location")) %>%
rename("species_richness" = "value")
#modA<-glmmTMB(species_richness ~ U_landscape+U_local+(1|plotid), data=df_species_rich, family=poisson) #poisson (no overdispersion)
#summary(modA)
#Anova(modA, type="3") #U_local is significant
#plot(simulateResiduals(modA))
#testDispersion(simulateResiduals(modA))
#plot(allEffects(modA))
dip_sp<-ggplot(df_species_rich, aes(x = U_landscape, y = species_richness, fill = U_local)) +
geom_boxplot(position=position_dodge(width=0.75)) +
geom_point(position=position_dodge(width=0.75), colour="gray3")+
scale_colour_manual(values = urb.col3) +
scale_fill_manual(values=urb.col3)+
theme_bw()+
#theme(legend.position = "none",
labs(x = "urbanisation level at local scale",
y = "Number of functional species per subplot",
title = "Species richness funct species DIPTERA")
#dip_sp
### graph resembles previous because those orders make up the majority of the species
## new way:
NMDS.vect <- scores(nmds_matrix, display = "species") %>%
as.data.frame()
NMDS.scores <- scores(nmds_matrix, display = "sites") %>%
as.data.frame()%>%
rownames_to_column("location") %>%
full_join(info, by="location")
stressplot(nmds_matrix)
# envfit() takes the output of metaMDS() and the species matrix you created
fit <- envfit(nmds_matrix, dip.matrix, perm = 999)
#extract p-values for each species
fit_pvals <- fit$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_spp <- fit %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
full_join(., fit_pvals, by = "func_species") %>%
filter(pvals <= 0.001)
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_local, shape = U_landscape), size = 3, alpha = 0.8) +
stat_ellipse(aes(color = U_local)) +
scale_colour_manual(name= "urbanisation level \nat local scale", values = urb.col3)+
scale_shape_manual(name= "urbanisation level \nat landscape scale", values=c(16, 17, 15))+
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.2, height=0.1)) +
clean_background
nmds_plot_new
###permanova diptera
library(pairwiseAdonis)
ST.dist<-vegdist(dip.matrix, distance="bray")
perman <- adonis2(ST.dist~ plotid + U_landscape * U_local,
data=info,
strata=info$plotid,
permutations=999)
perman
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = ST.dist ~ plotid + U_landscape * U_local, data = info, permutations = 999, strata = info$plotid)
## Df SumOfSqs R2 F Pr(>F)
## plotid 26 5.9884 0.39988 1.4151 0.142
## U_local 2 0.6420 0.04287 1.9721 0.018 *
## U_landscape:U_local 4 0.5326 0.03556 0.8180 0.779
## Residual 48 7.8127 0.52169
## Total 80 14.9757 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(perman)
| Df | SumOfSqs | R2 | F | Pr(>F) | |
|---|---|---|---|---|---|
| plotid | 26 | 5.9884451 | 0.3998777 | 1.4150825 | 0.142 |
| U_local | 2 | 0.6419672 | 0.0428673 | 1.9720772 | 0.018 |
| U_landscape:U_local | 4 | 0.5325962 | 0.0355640 | 0.8180486 | 0.779 |
| Residual | 48 | 7.8126829 | 0.5216910 | NA | NA |
| Total | 80 | 14.9756915 | 1.0000000 | NA | NA |
bd_local <- betadisper(ST.dist, info$U_local)
bd_local
##
## Homogeneity of multivariate dispersions
##
## Call: betadisper(d = ST.dist, group = info$U_local)
##
## No. of Positive Eigenvalues: 43
## No. of Negative Eigenvalues: 37
##
## Average distance to median:
## LOW MEDIUM HIGH
## 0.4067 0.3960 0.4131
##
## Eigenvalues for PCoA axes:
## (Showing 8 of 80 eigenvalues)
## PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
## 4.3459 1.4790 1.3238 1.1250 0.8922 0.7949 0.6316 0.5680
print(permutest(bd_local,pairwise=TRUE,permutations=999)) ## als het niet significant is, dan is het goed
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.00405 0.0020243 0.1404 999 0.868
## Residuals 78 1.12432 0.0144143
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.76200 0.852
## MEDIUM 0.75211 0.595
## HIGH 0.84450 0.58989
bd_land <- betadisper(ST.dist, info$U_landscape)
print(permutest(bd_land,pairwise=TRUE,permutations=999))
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.02294 0.011472 0.9652 999 0.385
## Residuals 78 0.92706 0.011885
##
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
## LOW MEDIUM HIGH
## LOW 0.35600 0.117
## MEDIUM 0.37465 0.736
## HIGH 0.11768 0.72640
no signs of “dispersion” => betadisper test is ok
info$interaction_comb <- interaction(info$U_landscape, info$U_local)
info$interaction_comb <- as.factor(info$interaction_comb)
pairwise.adonis(ST.dist, info$U_local) |> as_tibble() |> knitr::kable()
| pairs | Df | SumsOfSqs | F.Model | R2 | p.value | p.adjusted | sig |
|---|---|---|---|---|---|---|---|
| LOW vs MEDIUM | 1 | 0.3239383 | 1.7890322 | 0.0332602 | 0.058 | 0.174 | |
| LOW vs HIGH | 1 | 0.5164513 | 2.7520290 | 0.0502635 | 0.007 | 0.021 | . |
| MEDIUM vs HIGH | 1 | 0.1225613 | 0.6713256 | 0.0127456 | 0.749 | 1.000 |
pairwise.adonis(ST.dist, info$U_landscape)
## pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
## 1 HIGH vs MEDIUM 1 0.1015874 0.5156143 0.009818305 0.898 1.000
## 2 HIGH vs LOW 1 0.3340548 1.8474924 0.034309721 0.053 0.159
## 3 MEDIUM vs LOW 1 0.2306406 1.2737888 0.023910235 0.200 0.600
atypical sample P24SG
data.st.dip <- data.st %>%
dplyr::filter(location != "P24SG" & order == "HEMIPTERA") %>%
droplevels()
dip.matrix <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat)) %>%
column_to_rownames(var="location")
##dit is de community matrix
##sum(dip.matrix) gives now 5287
info <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat)
## dit is de info
nmds_matrix=metaMDS(dip.matrix, distance="bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2880867
## Run 1 stress 0.2848808
## ... New best solution
## ... Procrustes: rmse 0.06113589 max resid 0.2215116
## Run 2 stress 0.2902554
## Run 3 stress 0.2967525
## Run 4 stress 0.2825922
## ... New best solution
## ... Procrustes: rmse 0.06269567 max resid 0.3127576
## Run 5 stress 0.2802908
## ... New best solution
## ... Procrustes: rmse 0.05111893 max resid 0.2331626
## Run 6 stress 0.2837882
## Run 7 stress 0.2818612
## Run 8 stress 0.2890104
## Run 9 stress 0.2854956
## Run 10 stress 0.2974851
## Run 11 stress 0.2834486
## Run 12 stress 0.2917573
## Run 13 stress 0.2889322
## Run 14 stress 0.2921012
## Run 15 stress 0.2930148
## Run 16 stress 0.2817663
## Run 17 stress 0.2832978
## Run 18 stress 0.2819635
## Run 19 stress 0.2955729
## Run 20 stress 0.2804275
## ... Procrustes: rmse 0.04593515 max resid 0.2261187
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 4: no. of iterations >= maxit
## 16: stress ratio > sratmax
nmds_matrix
##
## Call:
## metaMDS(comm = dip.matrix, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dip.matrix))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2802908
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 5 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dip.matrix))'
sppr<- specnumber(dip.matrix)
sppr
## P01SG P01SY P01SR P02SG P02SY P02SR P03SG P03SY P03SR P04SG P04SY P04SR P05SG
## 6 5 11 4 9 4 5 7 7 4 9 6 13
## P05SY P05SR P06SG P06SY P06SR P07SG P07SY P07SR P08SG P08SY P08SR P09SG P09SY
## 11 7 8 6 4 9 9 7 8 7 4 9 5
## P09SR P10SG P10SY P10SR P11SG P11SY P11SR P12SG P12SY P12SR P13SG P13SY P13SR
## 4 14 7 9 12 5 5 8 11 8 11 3 10
## P14SG P14SY P14SR P15SG P15SY P15SR P16SG P16SY P16SR P17SG P17SY P17SR P18SG
## 10 7 6 5 12 5 8 9 12 5 5 8 10
## P18SY P18SR P19SG P19SY P19SR P20SG P20SY P20SR P21SG P21SY P21SR P22SG P22SY
## 7 5 9 7 6 5 8 6 9 5 15 7 5
## P22SR P23SG P23SY P23SR P24SY P24SR P25SG P25SY P25SR P26SG P26SY P26SR P27SG
## 13 8 6 5 2 5 5 8 3 9 6 5 4
## P27SY P27SR
## 4 3
df_species_rich <- sppr %>%
enframe() %>%
full_join(info, by = c("name"="location")) %>%
rename("species_richness" = "value")
modA<-glmmTMB(species_richness ~ U_landscape+U_local+(1|plotid), data=df_species_rich, family=poisson) #poisson
summary(modA)
## Family: poisson ( log )
## Formula: species_richness ~ U_landscape + U_local + (1 | plotid)
## Data: df_species_rich
##
## AIC BIC logLik deviance df.resid
## 392.1 406.4 -190.0 380.1 74
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plotid (Intercept) 4.72e-09 6.87e-05
## Number of obs: 80, groups: plotid, 27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.98105 0.09310 21.279 <2e-16 ***
## U_landscapeMEDIUM 0.09624 0.10461 0.920 0.358
## U_landscapeHIGH 0.15094 0.10222 1.477 0.140
## U_localMEDIUM -0.14086 0.10143 -1.389 0.165
## U_localHIGH -0.15173 0.10172 -1.492 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3") #U_local is significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: species_richness
## Chisq Df Pr(>Chisq)
## (Intercept) 452.7930 1 <2e-16 ***
## U_landscape 2.2090 2 0.3314
## U_local 2.8254 2 0.2435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0533, p-value = 0.656
## alternative hypothesis: two.sided
plot(allEffects(modA))
ggplot(df_species_rich, aes(x = U_landscape, y = species_richness, fill = U_local)) +
geom_boxplot(position=position_dodge(width=0.75)) +
geom_point(position=position_dodge(width=0.75), colour="gray3")+
scale_colour_manual(values = urb.col3) +
scale_fill_manual(values=urb.col3)+
theme_bw()+
#theme(legend.position = "none",
labs(x = "urbanisation level at local scale",
y = "Number of functional species per subplot",
title = "Species richness funct species HEMIPTERA")
### graph resembles previous because those orders make up the majority of the species
## new way:
NMDS.vect <- scores(nmds_matrix, display = "species") %>%
as.data.frame()
NMDS.scores <- scores(nmds_matrix, display = "sites") %>%
as.data.frame()%>%
rownames_to_column("location") %>%
full_join(info, by="location")
stressplot(nmds_matrix)
fig_nmds <- ggplot() +
geom_point(data = NMDS.scores, aes(x = NMDS1, y = NMDS2 , color = U_landscape)) +
scale_color_manual(values = urb.col3) +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
clean_background +
labs(x = "NMDS1",
y = "NMDS2",
title = "NMDS func groups")
fig_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_landscape, shape = U_landscape)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
stat_ellipse(linetype = 2, size = 1) +
clean_background +
labs(title = "NMDS")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
plot_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_local, shape = U_local)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
stat_ellipse(linetype = 2, size = 1) +
clean_background +
labs(title = "NMDS")
plot_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_local, shape =U_landscape)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
clean_background +
stat_ellipse(linetype =2 , size=1)+
labs(title = "NMDS")
plot_nmds
# envfit() takes the output of metaMDS() and the species matrix you created
fit <- envfit(nmds_matrix, dip.matrix, perm = 999)
#extract p-values for each species
fit_pvals <- fit$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_spp <- fit %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
full_join(., fit_pvals, by = "func_species") %>%
filter(pvals <= 0.001)
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_landscape, shape = U_landscape), size = 3, alpha = 0.8) +
stat_ellipse(aes(color = U_landscape)) +
scale_color_manual(values = urb.col3) +
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.3, height=0.1)) +
clean_background
nmds_plot_new
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_local, shape = U_landscape), size = 3, alpha = 0.8) +
scale_colour_manual(name= "Urbanisation \nat local scale", values= urb.col3)+
scale_shape_manual(name= "Urbanisation \nat local scale", values=c(16, 17, 15))+
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.2, height=0.1)) +
clean_background
nmds_plot_new
atypical sample P20SG
data.st.dip <- data.st %>%
dplyr::filter(location !="P20SG" & order == "HYMENOPTERA") %>%
droplevels()
dip.matrix <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat)) %>%
column_to_rownames(var="location")
##dit is de community matrix
##sum(dip.matrix) gives now 5287
info <- data.st.dip %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat)
## dit is de info
nmds_matrix=metaMDS(dip.matrix, distance="bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2456012
## Run 1 stress 0.259548
## Run 2 stress 0.2471578
## Run 3 stress 0.257863
## Run 4 stress 0.2569904
## Run 5 stress 0.2526869
## Run 6 stress 0.253721
## Run 7 stress 0.2455937
## ... New best solution
## ... Procrustes: rmse 0.003113757 max resid 0.02404952
## Run 8 stress 0.2659906
## Run 9 stress 0.2558501
## Run 10 stress 0.2577476
## Run 11 stress 0.2544416
## Run 12 stress 0.2506899
## Run 13 stress 0.246599
## Run 14 stress 0.2455937
## ... New best solution
## ... Procrustes: rmse 0.0001967132 max resid 0.0009612624
## ... Similar to previous best
## Run 15 stress 0.2566664
## Run 16 stress 0.2540408
## Run 17 stress 0.2524546
## Run 18 stress 0.2523057
## Run 19 stress 0.2503778
## Run 20 stress 0.2507785
## *** Best solution repeated 1 times
nmds_matrix
##
## Call:
## metaMDS(comm = dip.matrix, distance = "bray")
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(dip.matrix))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.2455937
## Stress type 1, weak ties
## Best solution was repeated 1 time in 20 tries
## The best solution was from try 14 (random start)
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(dip.matrix))'
sppr<- specnumber(dip.matrix)
sppr
## P01SG P01SY P01SR P02SG P02SY P02SR P03SG P03SY P03SR P04SG P04SY P04SR P05SG
## 2 4 2 7 3 5 2 5 4 4 7 4 6
## P05SY P05SR P06SG P06SY P06SR P07SG P07SY P07SR P08SG P08SY P08SR P09SG P09SY
## 4 3 10 3 3 6 9 5 8 2 7 10 6
## P09SR P10SG P10SY P10SR P11SG P11SY P11SR P12SG P12SY P12SR P13SG P13SY P13SR
## 5 10 3 4 8 5 7 8 8 3 5 3 6
## P14SG P14SY P14SR P15SG P15SY P15SR P16SG P16SY P16SR P17SG P17SY P17SR P18SG
## 8 6 8 6 9 3 4 6 3 7 8 7 5
## P18SY P18SR P19SG P19SY P19SR P20SY P20SR P21SG P21SY P21SR P22SG P22SY P22SR
## 3 8 3 4 4 3 9 4 4 4 1 4 5
## P23SG P23SY P23SR P24SG P24SY P24SR P25SG P25SY P25SR P26SG P26SY P26SR P27SG
## 7 6 7 3 3 3 5 1 5 3 3 6 6
## P27SY P27SR
## 5 9
df_species_rich <- sppr %>%
enframe() %>%
full_join(info, by = c("name"="location")) %>%
rename("species_richness" = "value")
modA<-glmmTMB(species_richness ~ U_landscape+U_local+(1|plotid), data=df_species_rich, family=nbinom2) #nbinom 2 to correct for overdispersion
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
summary(modA)
## Family: nbinom2 ( log )
## Formula: species_richness ~ U_landscape + U_local + (1 | plotid)
## Data: df_species_rich
##
## AIC BIC logLik deviance df.resid
## 361.8 378.5 -173.9 347.8 73
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plotid (Intercept) 0.00168 0.04099
## Number of obs: 80, groups: plotid, 27
##
## Dispersion parameter for nbinom2 family (): 4.29e+07
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.81972 0.10519 17.299 <2e-16 ***
## U_landscapeMEDIUM -0.10392 0.11938 -0.870 0.384
## U_landscapeHIGH -0.15392 0.12234 -1.258 0.208
## U_localMEDIUM -0.18834 0.12102 -1.556 0.120
## U_localHIGH -0.09803 0.11818 -0.830 0.407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(modA, type="3") #U_local is significant
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: species_richness
## Chisq Df Pr(>Chisq)
## (Intercept) 299.2621 1 <2e-16 ***
## U_landscape 1.6882 2 0.4299
## U_local 2.4331 2 0.2963
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(simulateResiduals(modA))
testDispersion(simulateResiduals(modA))
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.92833, p-value = 0.696
## alternative hypothesis: two.sided
plot(allEffects(modA))
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
## Warning in Effect.glmmTMB(predictors, mod, vcov. = vcov., ...): overriding
## variance function for effects/dev.resids: computed variances may be incorrect
ggplot(df_species_rich, aes(x = U_landscape, y = species_richness, fill = U_local)) +
geom_boxplot(position=position_dodge(width=0.75)) +
geom_point(position=position_dodge(width=0.75), colour="gray3")+
scale_colour_manual(values = urb.col3) +
scale_fill_manual(values=urb.col3)+
theme_bw()+
#theme(legend.position = "none",
labs(x = "urbanisation level at local scale",
y = "Number of functional species per subplot",
title = "Species richness funct species HYMENOPTERA")
### graph resembles previous because those orders make up the majority of the species
## new way:
NMDS.vect <- scores(nmds_matrix, display = "species") %>%
as.data.frame()
NMDS.scores <- scores(nmds_matrix, display = "sites") %>%
as.data.frame()%>%
rownames_to_column("location") %>%
full_join(info, by="location")
stressplot(nmds_matrix)
fig_nmds <- ggplot() +
geom_point(data = NMDS.scores, aes(x = NMDS1, y = NMDS2 , color = U_landscape)) +
scale_color_manual(values = urb.col3) +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
clean_background +
labs(x = "NMDS1",
y = "NMDS2",
title = "NMDS func groups")
fig_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_landscape, shape = U_landscape)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
stat_ellipse(linetype = 2, size = 1) +
clean_background +
labs(title = "NMDS")
plot_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_local, shape = U_local)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
stat_ellipse(linetype = 2, size = 1) +
clean_background +
labs(title = "NMDS")
plot_nmds
plot_nmds <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2, color = U_local, shape =U_landscape)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_manual(values = urb.col3) +
clean_background +
stat_ellipse(linetype =2 , size=1)+
labs(title = "NMDS")
plot_nmds
# envfit() takes the output of metaMDS() and the species matrix you created
fit <- envfit(nmds_matrix, dip.matrix, perm = 999)
#extract p-values for each species
fit_pvals <- fit$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_spp <- fit %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("func_species") %>%
full_join(., fit_pvals, by = "func_species") %>%
filter(pvals <= 0.001)
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_landscape, shape = U_landscape), size = 3, alpha = 0.8) +
stat_ellipse(aes(color = U_landscape)) +
scale_color_manual(values = urb.col3) +
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.3, height=0.1)) +
clean_background
nmds_plot_new
nmds_plot_new <- ggplot(NMDS.scores, aes(x = NMDS1, y = NMDS2)) +
coord_fixed() +
geom_point(aes(color = U_local, shape = U_landscape), size = 3, alpha = 0.8) +
scale_colour_manual(name= "Urbanisation \nat local scale", values= urb.col3)+
scale_shape_manual(name= "Urbanisation \nat local scale", values=c(16, 17, 15))+
geom_segment(data = fit_spp, aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm")),
col = "black") +
geom_text(data = fit_spp, aes(label = func_species), position = position_jitter(width = 0.2, height=0.1)) +
clean_background
nmds_plot_new
#biodiversity indices
data.st.corr <- data.st
#%>% filter(location!="P20SG" & location!="P24SG") %>%
# droplevels() ## dit verandert er niets aan dus heb ik niet uitgevoerd
func.matrix <- data.st.corr %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
subset(select=-c(plotid, U_landscape, U_local, urb_cat)) %>%
column_to_rownames(var="location")
##dit is de community matrix
info <- data.st.corr %>%
dplyr::group_by(location, plotid, U_landscape, U_local, urb_cat, func_species) %>%
dplyr::summarise(total_n_indiv=sum(n_indiv)) %>%
spread(key=func_species, value=total_n_indiv, fill=0) %>%
select(location, plotid, U_landscape, U_local, urb_cat)
## dit is de info
sppr<-specnumber(func.matrix) ## will give number of prey order within each sample
shannondiv <- diversity(func.matrix, index=c("shannon"))
#Pielou's evenness
Pielou <- shannondiv/log(sppr)
df_div <- data.frame(sppr, shannondiv, Pielou)
library(tibble)
df_div <- rownames_to_column(df_div, var="location")
df_div <- left_join(df_div, info, by="location")
df_div$U_landscape <- factor(df_div$U_landscape, levels = c("LOW", "MEDIUM", "HIGH"))
df_div$U_local <- factor(df_div$U_local, levels = c("LOW", "MEDIUM", "HIGH"))
sp_rich<-ggplot(df_div, aes(U_local, sppr, fill= U_local)) +
geom_boxplot()+
geom_jitter(width=0.25)+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
facet_grid(~U_landscape)+
ylab("Species Richness")+
xlab("")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1.1)))
sp_rich
g1 <- ggplot_gtable(ggplot_build(sp_rich))
stript <- which(grepl('strip-t', g1$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g1$grobs[[i]]$grobs[[1]]$childrenOrder))
g1$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g1)
sp_shannon<-ggplot(df_div, aes(U_local, shannondiv, fill= U_local)) +
geom_boxplot()+
geom_jitter(width=0.25)+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
facet_grid(~U_landscape)+
ylab("shannon diversity")+
xlab("")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1.1)))
sp_shannon
g_sh <- ggplot_gtable(ggplot_build(sp_shannon))
stript <- which(grepl('strip-t', g_sh$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g_sh$grobs[[i]]$grobs[[1]]$childrenOrder))
g_sh$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g_sh)
sp_pi<-ggplot(df_div, aes(U_local, Pielou, fill= U_local)) +
geom_boxplot()+
geom_jitter(width=0.25)+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
facet_grid(~U_landscape)+
ylab("Pielou Eveness")+
xlab("")+
theme_bw()+
theme(axis.text = element_text(size=rel(1)), axis.title=element_text(size=rel(1.1)))
sp_pi
g_pi <- ggplot_gtable(ggplot_build(sp_pi))
stript <- which(grepl('strip-t', g_pi$layout$name))
fills <- scales::alpha(urb.col3, 0.30)
k <- 1
for (i in stript) {
j <- which(grepl('rect', g_pi$grobs[[i]]$grobs[[1]]$childrenOrder))
g_pi$grobs[[i]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1
}
ggarrange(g_pi)
fig_div<-ggarrange(g1, g_sh, g_pi, nrow=3)
fig_div
annotate_figure(fig_div,bottom=text_grob("facets = Urbanisation level at landscape scale", vjust=0.2, size=15))
d1<-ggplot(df_div, aes(U_landscape, sppr, fill= U_local)) +
geom_boxplot(position=position_dodge(width=0.75))+
geom_jitter(position=position_dodge(width=0.75), alpha=0.6, size=1.5, color="gray3")+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
ylab("Species Richness")+
xlab("urbanisation level at landscape scale")+
theme_bw()+
theme(axis.text = element_text(size=rel(0.9)), axis.title=element_text(size=rel(1)))
d2<-ggplot(df_div, aes(U_landscape, shannondiv, fill= U_local)) +
geom_boxplot(position=position_dodge(width=0.75))+
geom_jitter(position=position_dodge(width=0.75), alpha=0.6, size=1.5, color="gray3")+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
ylab("shannon diversity")+
xlab("urbanisation level at landscape scale")+
theme_bw()+
theme(axis.text = element_text(size=rel(0.9)), axis.title=element_text(size=rel(1)))
d3<-ggplot(df_div, aes(U_landscape, Pielou, fill= U_local)) +
geom_boxplot(position=position_dodge(width=0.75))+
geom_jitter(position=position_dodge(width=0.75), alpha=0.6, size=1.5, color="gray3")+
scale_fill_manual(values = urb.col3, name="urbanisation level \nat local scale") +
ylab("Pielou Eveness")+
xlab("urbanisation level at landscape scale")+
theme_bw()+
theme(axis.text = element_text(size=rel(0.9)), axis.title=element_text(size=rel(1)))
fig_d <- ggarrange(d1, d2, d3, nrow=3)
annotate_figure(fig_d,top=text_grob("Effect of urbanisation on functional diversity arthropod communities", size=15, hjust=0.65))